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ABSTRACT 


A method for solving the viscous shock-layer equations for 
hypersonic flows over long slender bodies is presented. The governing 
equations are solved by employing a spatial-marching implicit finite- 
difference technique. The two first-order equations, continuity and 
normal momentum, are solved simultaneously as a coupled set. This 
method yields a simple and computationally efficient technique. 

Flows past hyperboloids and sphere cones with body half angles or 
five to 35 degrees are considered. The flow conditions included are 
from high Reynolds numbers at low altitudes to low Reynolds numbers at 
high altitudes. Detailed comparisons have been made with other 
predictions and experimental data for slender body flows. 

The results show that the coupling between the continuity and 
normal momentum equations is essential and adequate to obtain stable and 
accurate solutions past long slender bodies. Both the Cebeci-Smith and 
Baldwin-Lomax turbulence models are found to be adequate for application 
to long slender bodies. Using the corrected slip models, the viscous 
shock-layer predictions compare quite favorably with experimental data. 
Under chemical nonequilibrium conditions, the surface catalytic effects 
can significantly influence the surface heat transfer. 
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Chapter 1 
INTRODUCTION 

1.1 Physical Features of the High Mach Number Flow 
A renewed interest in hypersonic aerothermodynamics has been 
motivated by new vehicle concepts such as the national aero-space plane 
* 

(NASP) [1] , and the transatmospheric vehicle (TAV) [2]. The term 
"hypersonic" implies that the flight velocity is much greater than the 
ambient speed of sound. An approximate classification of this flow 
regime is where the freestream Mach number is greater than five. These 
vehicles will encounter a variety of flow conditions which include 
atmospheric flight and the transitional flow regimes. For example [ 3 ], 
the TAV will take off from the Earth’s surface, and enter a low Earth 
orbit. The vehicle will carry out a global mission either inside or 
outside the atmosphere, and eventually land back on Earth under its own 
power. The range of flow conditions includes low altitude, high density 
flow, to high altitude, low density flow. These conditions include the 
continuum flow regime where the no-slip assumption is made. In addition 
the transitional flow regime, where slip effects are important, must be 
considered. 

*The numbers in brackets indicate references. 
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There are two main effects associated with hypersonic flows [4,5]: 

(1) the fluid dynamic effects arising from the high velocity gas, and 

(2) the real gas effects due to the high temperature gas. For 
hypersonic speeds, the shock layer, which is defined as the distance 
between the shock and body, is small. The boundary layer thickness 
grows more rapidly because kinetic energy dissipation within the 
boundary layer, which yields a high gas temperature, results in an 
increase in gas viscosity and a decrease in density. Along with the 
thin shock layer, the thick boundary layer creates an important 
disturbance in the outer flow that gives rise to the viscous interaction 
phenomenon which controls the surface pressure distribution over the 
body. Moreover, for a blunt body, the shock wave is curved, leading to 
large entropy gradients in the shock layer. 

Compression of the gas forward of the vehicle and heat generation 
due to viscous dissipation lead to elevate gas temperatures in the shock 
layer. Additionally, the gas will promote chemical reaction in both the 
boundary layer and the shock layer. As a result, the specific heat per 
unit mass is increased considerably, and the specific heat ratios will 
no longer be constant. The gas will not behave as a calorically perfect 
gas. Moreover, if the shock layer temperature is sufficiently elevated, 
radiation effects will become important, giving rise to a radiative flux 
to the surface. The resulting heat transfer to the surface of a 
hypersonic vehicle will dominate the design criteria of the vehicle. 

For hypersonic speeds, Blottner [6] showed that the shock layer 
flow is in chemical equilibrium and has a definite boundary layer region 
for low altitude conditions. Also, the flow can become turbulent at 
these conditions [7]. However, the gas may not reach the equilibrium 



3 


state for higher altitude conditions. Moreover, the boundary layer 
cannot be identified because it merges with the shock wave. 

1.2 Numerical Methods for Hypersonic Flow 
In general, there are two methods to analyse hypersonic flow — 
experimental and theoretical . Flows in chemical equilibrium can be 
simulated with small-scale laboratory experiments with corrections for 
"real gas" effects. However, the chemical nonequilibrium flow around a 
hypersonic vehicle operating in the upper atmosphere cannot be simulated 
because it requires simultaneous reproduction of air density, flight 
velocity and vehicle scale. In absence of a full-scale flight 
experiment in which the thermodynamic environment is fully duplicated, 
an adequate design capability for hypersonic vehicles relies on 
theoretical predictions. 

The great entropy gradients and the thick boundary layer in a 
hypersonic flow make the classical isentropic irrotational approach and 
the conventional first-order boundary layer equations inadequate to 
predict the flowfield. Second-order boundary layer effects and 
vorticity interaction should be considered in the flow. Three current 
numerical approaches have been adopted for analysing these problems. 

They are the solution of either the second-order boundary layer 
equations, the full Navier-Stokes equations, or the viscous shock-layer 
equations. 

The simplest of these approaches is to employ the second-order 
boundary layer equations [8], with matching of the first- and second- 
order boundary layer and inviscid solutions at the boundary layer edge. 
Although this approach has been found quite attractive for short slender 
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bodies [9], several computational difficulties become obvious for long 
slender bodies. First, the computing time required is excessive because 
one must compute the inviscid flow, the first-order boundary layer flow, 
the flow due to displacement thickness, and then the second order 
boundary layer flow. A second difficulty arises from strong vorticity 
interaction which may occur far downstream due to the entropy layer on 
blunt long slender bodies. This entropy layer causes difficulty in the 
matchup procedure between the viscous and inviscid regions because there 
is a question as to what the proper edge conditions should be for the 
boundary layer. This approach does not properly take into account the 
swallowing of the strong entropy layer by the boundary layer [5] . The 
design geometries of hypersonic vehicles are slender long bodies with 
blunt noses in order to reduce heat transfer rate at the stagnation 
region and reduce drag force on the body. The second-order boundary 
layer equations are not desirable for this problem. 

The second approach employs the steady Navier-Stokes equations 
[10] and their time- dependent forms [11]. This approach successfully 
provides the solution for the stagnation region of short bodies. 

However, the complexity of the solution procedure due to the elliptic 
nature of the equations requires excessive computing time and computer 
storage, which currently limits their applications to short bodies. 

Because of the difficulties encountered by the above two 
approaches, attention has turned toward the third approach, the viscous 
shock-layer equations. This set of equations which was developed by 
Davis [12] is obtained from the full Navier-Stokes equations by keeping 
terms up to second-order in the inverse square root of the Reynolds 
number in both viscous and inviscid regions. It is uniformly valid to 
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second-order throughout the entire shock layer, hence, the viscous- 
inviscid interactions and strong vorticity interactions are accounted 
for in a straight-forward manner. Moreover, this set of equations is of 
a hyperbolic-parabolic nature and, therefore, can be solved by using a 
marching procedure similar to methods employed in boundary layer theory. 
As a consequence, they can be solved for a hypersonic flow on a slender 
long body without excessive computer time and storage requirements. 
Moreover, this set of equations can be used to compute the viscous flow 
in the subsonic blunt nose region. This is desirable for long bodies, 
especially for analysing problems with chemical reactions. 

The full viscous shock layer solution of Davis [12] was obtained 
through an iterative relaxation process from the thin shock layer 
solution. This approach encountered difficulties for the flow far 
downstream, especially for the slender body. Werle et al. [13] 
developed an Alternating Direction Implicit (ADI) Technique with an 
artificial time coordinate to relax the shock shape from an initial 
guess. Even with large relaxation factors, the instabilities were still 
encountered whenever the inviscid region encompasses a significant 
portion of the total shock layer thickness. The relaxations of the 
shock shape in Davis [12] and Werle et al. [13] are essential due to the 
slightly elliptic nature of the equations in the streamwise direction. 
These instabilities do not come from the shock shape relaxing technique. 

From the hypersonic small disturbance theory [14], it is shown 
that the continuity, normal momentum, and energy equations become 
uncoupled from the tangential momentum equation in the inviscid region. 
In other words, the solution of the continuity, normal momentum and 
energy equations will not depend strongly on the solution of the 
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tangential momentum equation. The numerical methods of Davis [12] and 
Werle et al. [13] solved the governing equations separately which are 
known as the cascading method. In this method, the solution of the 
tangential momentum equation drives the solution of the continuity and 
normal momentum equations. This method becomes improper for the flow 
far downstream, especially for a slender body on which the shock layer 
thickness is very thick and the inviscid region encompasses a large 
portion of the shock layer. An alternative method of solution was 
suggested by Werle et al. [13]. The more adequate method is to solve 
the equations simultaneously. 

A fully coupled system of all the equations is a desirable scheme. 
Hosny et al. [15] solved the four governing equations, namely, the 
continuity, tangential momentum, normal momentum and energy equations, 
simultaneously as a coupled set and local iterations were made to solve 
for the shock stand-off distance. Gorden and Davis [16] added an 
equation for the shock stand-off distance into the coupled set to 
eliminate the need for local iterations. This technique is quite 
appealing for perfect gas applications. But, it requires inversion of 
large matrices and hence the storage and computing requirements are 
quite large. Also, the system of equations will become very complicated 
if chemical reactions are included. Therefore, this approach is not 
desirable for long bodies. 

The two second-order equations, tangential momentum and energy, 
are parabolic, and there are few problems in finding the solutions to 
these equations. The greatest difficulty exists in solving the two 
first-order equations, continuity and normal momentum [13]. Moreover, 
from the hypersonic small disturbance theory, the solution of the 
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tangential momentum equation becomes uncoupled with the other equations 
far downstream for a slender body. A more desirable approach is to 
solve the two first-order equations simultaneously as a coupled set 
rather than solve all four equations as a fully coupled set. This 
approach is quite attractive for slender body problems, especially with 
real gas effects. The instabilities will be eliminated, and the storage 
requirements and computing time may not increase excessively. 

Waskiewicz and Lewis [17] coupled the two first-order equations and 
reported good improvement in the solution obtained for slender (7 degree 
and 10 degree) but short bodies with solutions up to 20 nose radii or 
less. The effects of this technique on the flow field over a slender 
long body far downstream should be investigated. 

Most of the work with viscous shock-layer equations in the past 
has considered either short slender or wide angle bodies. However, most 
of the future vehicles will be long slender blunt bodies. The 
calculation of hypersonic viscous flows past long slender axisymmetric 
blunt bodies is of prime interest to the designer of such aerospace 
vehicles. 

A variety of flow conditions are encountered during the 
transatmospheric flight of these vehicles. The range is from low 
Reynolds numbers at high altitudes to high Reynolds numbers at low 
altitudes. At low altitude, the hypersonic flow over a slender body 
usually becomes turbulence. Direct numerical solution of turbulent flow 
cannot be obtained at the present. The prediction of turbulent effects 
depends on modeling the fluctuation terms. The algebraic eddy viscosity 
models are more appealing than the other models because less computer 
storage and computer time are required. The general algebraic 
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turbulence model which is implemented with the viscous shock layer 
method is the Cebeci-Smith model [18,193. Due to the difficulty of 
determining the boundary layer edge in a hypersonic flow over a long 
slender body, an alternative model, the Baldwin-Lomax model [20], is 
more likely to be used. This model has been implemented in the Navier- 
Stokes equations and the parabolic Navier-Stokes equations [21,22], but 
not in the viscous shock-layer equations. 

At high altitude, the "low density effects" become important where 
they can significantly influence the lift, drag, moments and aerodynamic 
heating of a hypersonic vehicle. However, not much attention has been 
given to the problems encountered with low-density aerothermodynamics. 

At highly rarefied gas flow conditions, the continuum approach is no 
longer valid. But at slightly rarefied gas flow conditions, the 
continuum approach can be extended to this flow regime if slip effects 
are properly accounted for [12,233. 

At high temperature, the perfect gas assumption becomes invalid. 
The gases in the flow become chemically reacting, especially, at high 
altitude. While flows with chemical equilibrium in the shock layer have 
been studied intensively [2M,25], there are only a limited number of 
analyses on flows with chemical nonequilibrium [6,26]. The effects of 
finite-rate chemical reactions are not fully understood yet. 

The main objective of this study, therefore, is to investigate the 
solution for long slender blunt body; i.e., the continuity and normal 
momentum equations are solved simultaneously as a coupled set. The 
basic theoretical formulations used in this study are provided in Chap. 
2. The suitability of the Cebeci-Smith and Baldwin-Lomax models for the 
hypersonic flow past long slender bodies at low altitudes is examined. 
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examined. It is presented in Chap. 3. The low-density effects using a 
perfect gas model are presented in Chap. 4. The investigation of a 
chemical nonequilibrium flow is presented in Chap. 5. 
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Chapter 2 

BASIC THEORETICAL FORMULATION 

The viscous shock-layer equations for a perfect gas [12] and for a 
multicomponent gas mixture [24] are presented in this chapter. The 
physical model and coordinate system for a body are shown in Fig. 2.1. 
The flow in the shock layer is assumed to be axisymmetic, steady, 
viscous and compressible. The shock-layer gas is assumed in local 
thermodynamic equilibrium. 

2.1 Governing Equations For a Perfect Gas 
The viscous shock— layer equations are obtained from compressible 
Navier-Stokes equations which are written in the body-oriented 
coordinate system shown in Fig. 2.1. These equations are 
nondimensionalized with variables which are of order one in the region 
near the body surface (boundary layer) for large Reynolds number. The 
same set of equations are then written in variables which are of order 
one in the essentially inviscid region outside the boundary layer. 

Terms in each set of equations up to second-order in the inverse square 
root of a Reynolds number are kept. A comparison of the two sets of 
equations is then made and one set of equations is obtained which is 
valid to second order in both the outer and inner regions. A solution 
to this set of equations is thus uniformly valid to second-order in the 
entire shock layer for arbitrary Y. Anderson and Moss [19] used the 
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eddy viscosity approximation to replace the Reynolds stresses and 
turbulent heat flux to find the solution of the turbulent flow. The 
viscous shock-layer equations for a perfect gas with turbulent flow are 


[19] 

Continuity: 


— [(r+ncose) J pu] + |-[( 1+nic) (r+ncos0)^pv] - 0 
3s <* n 


( 2 . 1 ) 


3 -momentum: 


, u 3u 3u uvK y + _L lE _ e 2 fi_r u ( 1+e + )iH 
1 1+n< 3s V 3n 1+mr 1+nic 3s l 3n L 3n 


_ uu* -i + ( 2k + _i£°®i_)[(i+ e + )— 
1+n*r 1+ntc r+ncose' L an 


yux 
3n 1 +n< 


]} 


( 2 . 2 ) 


n-momentum: 


0 (_ ILiv + v 3v _ JijC) + |E . o 
p M+nic 3s 3n 1+nx 3n 


(2.3) 


Energy : 


0 ( 1L M ♦ v M) - v iE + £H_X! £ 

p ' 1 +n»c 3s V 3n' 3n 1+nic 


>3 r M/ , ♦ P r . 3H . . , k Jcose . 

3n^Pr +e Pr,t 3n ^ 1+nic r+ncose^ 


x [— (1+e + -^— ) + <t>3 ) 

x L Pr^ Pr ,t ; 3n v 


(2.4) 


State: 


( 7 - 1 ) 


pT 


P 


Y 


(2.5) 
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where 


♦-P 

The molecular viscosity is given by the Sutherland's law as 

r 3/2 


9u n _ yu tc 
I +nic 


y - [(1+C)/(T+C)] T- 


where 


* ox 

C - C /(Y-1 ) M T 


and C is 110.33 K for air. 

The preceding equations are written in nondimensional forms 
the nondimensional variables are defined as 


* * 
3 ■ s /r n 


# * 

u - u /U 


# # 

P - P /p U 


* * 

n - n /R n 


* * 
r - r /R n 


# * 

v - v /U 


* * 

p - p /p„ 


# * 

T - T /T _ 
ref 


* x x 

K ■ M /U <T ref > 


^ ft 
" " " B N 


X X 

h - h /h 


X X 

C - C /C 
P P P»‘ 


Also, the dimensionless parameters that appear in Eqs. (2.1) to (1 
are defined as 


xx x 

Pr - y C /k 
P 


XX X 

Pr-.t - „ T C p /K T 


r “ (T ref',1/2 
e “ L * * » J 


* * * 
P® 


and 


+ xx 

e - y T /y 


( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 

where 

(2.9) 

. 6 ) 

( 2 . 10 ) 
( 2.11 ) 

( 2 . 12 ) 


(2.13) 
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To simplify the numerical computations, a coordinate 
transformation is applied to the viscous shock-layer equations. This 
transformation is accomplished by normalizing the normal coordinate with 
respect to the local shock standoff distance. Consequently , a constant 
number of finite-difference grid points between the body and shock are 
used. Also, the need for interpolation to determine shock shape and the 
addition of grid points in the normal direction is eliminated as the 


computation moves downstream. 

The transformed variables are 
s - 5 


n 


- n/n 


sh 


(2. 1M) 

( 2 . 15 ) 


The transformations relating the differential quantities are 


- « n* 3 

1 _ _ 3 _ ( 2 . 16 ) 

3s 35 n 3h 11 3 tT 


3 _ 

3n 


1 3 



( 2 . 17 ) 


and 

3 m _ 1 _ 3 ( 2 . 18 ) 

3n 2 n 2 3n 2 
sh 

where 

n . . ( 2 . 19 ) 

sh d£ 

In general, a variable grid spacing is used in the n~direction so 
that the grid spacing can be made small in the region of large 
gradients. Since the spacing is not uniform in the n direction, it is 
convenient to apply a transformation to the n coordinate so that the 
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governing equations can be solved on a uniformly spaced grid in the 
computational plane. This transformation is achieved by [26] 


n - g(- — ) - g(n) 
sh 

The stretching function g(n) is given by 

i _ r~ . (1-a) a _ r 6 — n( 2a + 1 ) + 1 1 t 

8 " C ( T^(2S*1)-1 ,] 


( 2 . 20 ) 


( 2.21 ) 


The first and second derivatives of Eq. (2.21) are expressed as 

dg _ ( 1 -a) (2a+1 ) f 1 1 . 

dn " £ n (F + 1 j lB“n(2a+1 ) + 1 ] [F+n(2a+1 )-1]' 

FT 


dfg _ (1-g)(2a+1 ) 2 { 1 1 

dn 2 4n(|rf) [e-n(2a+1)+1] 2 [?+n(2a+l )-1 ] 2 


(2.23) 


Equation (2.21) permits the mesh to be refined either near the body only 
(o - 0) or equally near both the body and bow 3hock (a - 1/2). The 
parameter g controls the amount of stretching. The coordinate n can be 
obtained from inverting Eq. (2.21) and is expressed as 


1~n~a 

,B+K 1_a 
1 - MPT' 

' • (STTT [1 - 6 <— 


~ 1 


-}] 


(2.24) 


( FT ) 


i-s 


+ i 


This transformation keeps the body at n * 0 and the shock at p ■ 1 with 
uniform mesh in the computational coordinate n. 

After the governing equations, Eqs. (2.1) to (2.6), are written in 


the transformed n coordinates, the second-order partial differential 
equations are expressed as [19] 
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9 2 w d 2 g/dn 2 + o 1 (dg/dn) 3W <* 2 

3n 2 (dg/dn) 2 ^ (dg/dn) 2 


a 3 a U 3W 

+ _ p + _ p ar 

(dg/dn) (dg/dn) 


The quantity W represents u for the s-momentum equation and H for 
energy equation. The coefficients to are written as. 


s-momentum, W - u: 


*1 + 

1 y(1+e ) 


3y(1+e ) dg + 


n . ic(1+2e ) 
sh 


jn a h c °3 0 


3n d " 0^3 h O(1* e + ) r + nn 3h CO30 


n ;h n sh npu 


n . pv 
sh K 


e 2 y(1+e + )(1+nn 3h K) e 2 y(1+e + ) 


n . k 
sh 


y(1+e )0 + nn sh ic) 


3y(1+e ) dg 

3n dn 


2 

n . k 
sh 


(1+e )(1+nn ah ic) 


n . icpv 

r < + JC038 1 Sh 

h+nn^K r + nn 3h cose e 2 (1+nn 3h Oy(1+e + ) 


n 


sh 


e 2 (1+nn ah K)y(1 + e ) ^ 


_ rl£ _ n s hT1 d S l£l 
+ . l 3E n gh dn 3n J 


2 . 25 ) 


the 


( 2 . 26 ) 


( 2 . 27 ) 


( 2 . 28 ) 
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“4 " 


n 3h pU 


e^(1+nn g ^<)u(1 + e ) 


Energy, W - H: 


a, - 


•jj— ( 1 +e 

Pr 


♦ Pr . 
Pr, t 


i_ [P_ (1+e + _P£_)] + 

dn 3n l ?r K e Pr,t ;j 


Kn ah 


t j^sh 0039 

r+nn .cose 
sh 


n'.n . npu 
sh sh K 


e2 pr ( 1 +e+ Pr7t ) ( 1 + ^ n sh ic) 


n . p v 
sh K 

e 2 i^(1 + e 


Pr' 


+ Pr . 
Pr,t 


a 


2 


- 0 


n 


a- - 


sh rdg |± + n f — J£ + _4£os0 j 

+ Pr . L dn 3n shM + nn . < r+nn . cos J Y 


3 M (1 + Pr . L dr? 3n shM + nn < r+nn cos 

Pr^ Pr f t' 


~ n . KVpu 

t V dg 3p sh 1 

2 — ~ ”2 — " J 

e dn 3n e ^ 1+nn sh K ^ 


n sh pu 


e2 p? (1+e+ Pr7t )(1+ ^ n sh ,c) 


( 2 . 29 ) 


( 2 . 30 ) 


( 2 . 31 ) 


( 2 . 32 ) 


( 2 . 33 ) 
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where 


H- [Pr - 1 
Pr L 


+ e 


Pr , p _ t _ n l _H_ 3u _ jiu 
pF7P Pr ’ i;j n sh dfi3n i + f& 


‘^sh* 


(2.34) 


The remaining first-order equations are written as 
Global Continuity: 

1fC“ S h (r *" n Sh 00S9)J,,Ul * W^ f<r *" n 3h 0039> 
x [( 1+nn gh ic)pv - n^ h npu]} - 0 


n-momentum: 


pu 3v _ n sh npU 3v dg + _gv 3v dg 
(1+ ^ n sh K) n sh (1+75h sh K) 3n d?T n sh 3,1 


pu 2 < + _L lE - o 

(l + TTn 3h O n sh 3n dfT 


( 2 . 36 ) 


State: 


p - pT (Y-1)/Y 


(2.37) 


2.2 Governing Equations for a Multicomponent Mixture 
The conservation equations that describe a reacting multicomponent 
gas mixture can be found in the literature [27,28]. The viscous shock- 
layer equations for a nonequilibrium multicomponent gas mixture are 
obtained from the conservation equations employing the same procedure as 
for the perfect gas [24]. The nondimensional forms of global 
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continuity, s-momentum, and n-momentum are written in the same form as 

Eqs. (2.1) to (2.3) except the eddy viscosity, e + , is set to zero. 

For a chemical nonequilibrium flow, the energy equation is 
formulated in terms of the temperature instead of total enthalpy. In 
addition to these conservation equations, the species continuity 
equation and equation of state are needed to complete the set of 
equations. The nondimensional forms of these equations are expressed as 
[24] 


Energy : 

PC ( 


II + v !I) . (_J±_ l£ + v l£) . 

p'1+nic 3s 3n v 1+nic 3s 3n' 


N 


. 2 r 3 # 3T» # < jcos0 . 3T r T n 3T j. 

^3n ^Sn 1+nic r+ncose K 3n ^ J i C p,i3n 


N 


,3u 


<u 

1+n< 


) 2 ] - l 


i-1 


h if 

n i i 


( 2 . 38 ) 


Species continuity: 


u 3C i 3C i 2 

p( 7niic 5T * v ■ 4 1 - Z ~J 

(1+me)(r+ncos0) J 


x t~ [ ( 1+n<) (r+ncos0 )^ J^]} 


(2.39) 


State: 


P “ PT <=sV-> 
M C 


(2.40) 
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The nondimensional variables are as defined by Eq. (2.9) and 


# # * 

* ■ * /u ref C p.* 


^ A fc ^ 

4 1 ■ *iV p - u - 


(2.41 ) 


J i ' J iV u *ef 


Using the same coordinate transformations as given for a perfect 
gas, Eqs. (2.38) and (2.39) are written in the form of Eq. (2.25). The 
W in this equation represents T for the energy equation and for the 

species continuity equation. The coefficients a 1 to a^, in this case, 

are given by 
Energy, W - T : 


1 3< dg ( 


n .< 
sh 


J n 3 h °° 39 


„ h N 3 

sh 


y j c 

“l ‘ k dir ’ 1 * I ' n 3h l( ” r+flh 3h c03e ic ^ 1 P,i 


n 3h n 3h 1 ’ C D OT 
2 — 

e ic(l+nn ah ic) 


n . pC v 
sh K p 

2 

e ic 


(2.42) 


(2.43) 


n =v,H 


sh_ ^_1_ 3u dg _ 


n . 3n dn 
sh 


, n 2 n H s . 

— — ) 2 l h,w 

1+nn sh K e 2 < i-1 


i i 


un 


sh 


9p ri 


sh 


u 




v dg 3p 


e^(Unn s h<) 35 e < 1+nVe V n gh dn 3n 


(2.44) 
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2 _ 
n . pC u 

. . - 3h P 

*H 2 — 

e K(1+nn gh ic) 


(2.45) 


Species continuity, W - C^: 


. _L !!h f! . V K , n 3h cose 

l 1 " PL. 3n dn + 1+TTn . k + r+TTn . cose^ 

i sn sn 


n . n' npu n . pv 

sh sh K sh K 


e 2 PL 1 (1+nn 3h K) e 2 PL i 


(2.46) 


a 2 “ ° 


(2.47) 


a 


3 


1 3PM i dg PM i 
PL^ 3n dfF + PL^ 


n . k 

( + 

^sh* 


2 

n .cose w.n . 
3h . i sh 

' + ? 

E PL, 


r+fjh .cose 
sn 


(2.48) 


n . pu 
stv 


e PL^i+nn k) 


where 


PL, 


lb n 


PM, - 


N s 3C dg 

U V AK k 

PF J, 1 
k-1 

k^i 


ik 


(2.49) 


(2.50) 


3n dn 


(2.51 ) 
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t Le i 


i-k 


Ab 


ik 


(2.52) 


M. 


Le,. - [ 4 Le ’ ik + t 1 
M 


w * N 
M, 3 


4) i Le -ij c j] i,<k 

M., j-1 


j^i 


In Eq. (2.52), Le,. are the multicomponent Lewis numbers and M is 


ik 

molecular weight which is given by 
“* 1 


M 


N C. 
s i 

y * 

L M 
i-1 i 


(2.53) 


The mass flux, due to concentration gradients can be written as [29] 


v N s 3C k dg 

Jl ' J, 4bi,< ^ 


The equation 
p - pT 


of state is given by 
* “ # * 


R 


/ M 


C 

P, 


00 


(2.54) 


(2.55) 


2.3 Boundary Conditions 

At low altitudes, slip effects are not important. No-slip and no- 
temperature-jump boundary conditions are used on the body surface. The 
wall temperature and enthalpy are specified as constant. The boundary 


conditions at the shock are calculated by using the Rankine-Hugoniot 
relations [12,24]. 
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At high altitudes, the continuum flow assumption breaks down in 
the region next to the wall. The no-slip and no-temperature jump 
boundary conditions are no longer valid. As such, the slip and 
temperature jump boundary conditions should be used. 

2.3.1 Boundary Conditions on the Wall 

Shidlovsky [30] has shown that at the body surface the velocity 
slip and temperature jump conditions are of the same order as the 
Knudsen number. The Knudsen number is defined as the ratio of the 
molecular mean free path in the gas to a characteristic dimension of the 
flowfield. The no-slip boundary conditions (which correspond to 
continuum conditions) are obtained when Kn •* 0. However, when the flow 
density decreases, the mean free path becomes long compared to the 
characteristic length in a region next to the wall. This region is 
called the Knudsen layer. Under this condition, the slip conditions 
should be used. 

The slip conditions are assumed to exist across the Knudsen layer. 
The net fluxes of momentum and energy at the outer edge of the Knudsen 
layer are equated to the difference between the incident and reflected 
fluxes at the wall. These fluxes are assumed to be constant across the 
Knudsen layer and are obtained from the moments of the distribution 
functions. The slip conditions are then obtained from the balance 
equations of these fluxes and are given by [31] 

Velocity slip: 

nr 2-9 2 **s , 3u _ _jcu_. 

u s “ IT £ / p p ^3n 1+ntc's 
V s s 


( 2 . 56 ) 
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Pressure slip: 


, „ - 2 

4 Y 2-6 e y g 

p = ' p - * Ijg ' r ^ ^ 


p 3T 
— (tt-) 
P S 3n 3 


(2.57) 


Temperature slip: 


1 IT Y 2-0 e 2 y g 3 T 

■ T w * I Js W't’pF SX‘35’ 


(2.58) 

'S r 3 3 

In the derivation of the relation for the temperature slip, it is 
assumed that the internal energy is frozen during the reflection from 
the wall. The parameter 1) is the accommodation coefficient and its 

value is taken to be unity in this study. 

In the computational plane, Eqs. (2.56) to (2.58) are expressed as 

[62] 

Velocity slip: 




2 

, 2-0 s e y s r 1 dg 3u _ 
L - Azr 




<u 


r3 


P n 3h dfT3n (1 ^ n sh K) s 
S K S 


(2.59) 


Pressure slip: 


2 

4 , Y w2-0. e y s £s ,J_ dg 3T\ 


Pg “ P + 


l p 3 "sh d?T 3n 


( 2 . 60 ) 


Temperature slip: 


„ _ 1 nr r Y i f 2-0-\ ef y 3 r 1 dg 3T^ 

T s “ T w + 2 J 2 Y-1 0 ^ Pr r-r- n _ h dn 3n J 


( 2 . 61 ) 


P s p s 3h 


2 . 3.2 Boundary Conditions at the Shock 

The boundary conditions at the shock are the modified Rankine- 


Hugoniot relations developed by Cheng [32]. The shock equations are 
obtained from one-dimensional Navier-Stokes equations which are written 



25 


in the shock-oriented coordinate system. The tangential derivatives are 
neglected as compared to the normal derivatives. Then, the shock 
equations in dimensional forms are expressed as 
Continuity : 

i-j(pV) - 0 (2.62) 

30 


Tangential-momentum: 

* * 30 * 3 , * 30 * 

P * — “ — (P —*> 
30 30 30 


Normal-momentum: 

* # 

# # 


„ 37 A 3p_ 3 r lo * 37 

P 7 — 5f + — nr - — »[(2p + X ) j] 
30 30 30 30 


Energy: 


* * * 3 T* * 3 p* 3 , * 3 T* v _ 

p C 7 — =■ - 7 -*~r “ — T (< — *) + 


30 


30 30 


30 


* o *2 

* * 37 ^ * 30 . 

(2p + X )( — y) + p ( — y) 


30 


30 


(2.63) 


(2.6*0 


(2.65) 


Integrating Eqs. (2.62) to (2.65) from just behind the shock wave 
to freestream, the shock conditions are obtained. The nondimensional 
forms are given by 
Continuity: 

p 7 » -sina (2.66) 

K sh sh 

Tangential-momentum: 

‘H'sh * “sh 31 "” ' l2 ' 67) 

Normal-momentum: 

p . - p + sina (sina + 7 . ) 

*sh ® sh 


( 2 . 68 ) 
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Energy: 


2 

t 


(-H 

v Pr 


3T , 

an' 


sh 


sing , JJY_ 
? *■ 

(Y+1) 


T sing - 3l | - ° (Q . - cosa)‘ 
3 h 2 sh 


sin 2 g + [(-^-) 
Y-1 


4(Y~1 ) -. 1_ 

(Y+1) 2 M 2 
00 


2 4 2 

(Y+1 ) M sin a 


(2.69) 


Since velocity components tangent and normal to the shock are not 
the same as those tangent and normal to the body surface, 
transformations are needed to relate these quantities. The 
transformations are given by 

U . - u . cos(o-9) + v . sin(a-e) (2.70) 

sh sh sh 

V . - - u . sin(a-e) + v . cos(g-0) (2.71) 

sh sh sh 

Also, the transformations between the body-oriented and shock-oriented 
coordinates are given by 

3 « s cos(a-0) + n sin(a-0) (2.72) 

h - - s sin(a-0) + n cos(a-0) (2.73) 

The derivative with respect to fl is related to (s,n) coordinate as 

■|= « cos(a-0) - sin(a-0) -r— (2. 7^) 

drl of! oS 

Consequently, the shock conditions in the computational plane are 

expressed as 

Continuity: 


^sh ^sh 


-sing 


(2.75) 
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Tangential-momentum: 


dn 


u 3h {[c°3 ( a- 9) + n sin(o-e) ] — ^ 


d£ 


n 3h dn 3n 


A 


(+)|| sin(a-0)} gh + Q gh sina - sina cosa 


Normal-momentum: 


p sh “ p « + 3inct ( sin « + 7 3h ) 


Energy: 


e2( Pr ) sh U CO3 ( a “ 0 ) + *1 IT sin(a-e)] || 


d£ 


"sh dP 3n 


- || 3in(a-e)} gh + T gh sina - 5i2“ (u gh - cosa) 2 


sing | _MY „,„2_ ^ 2 N 4(Y-1) n 1 


2 1 2 
(Y+1) 2 


sin 2 a + [(-j-) - 


] -L_ 

(Y+1 ) 2 M 2 


2 U 2 ‘ 
(Y+1 ) M sin a 


Equation of state: 


p . - Yp ./( Y— 1 )T . 
K sh *sh sh 


( 2 . 76 ) 


(2.77) 


( 2 . 78 ) 


(2.79) 


From Eqs. (2.76) and (2.78), it is noted that the first term in these 
two equations can be neglected at low altitudes where e is very small. 
On the other hand, this term becomes important when e is large. 



28 


Chapter 3 

HIGH REYNOLDS NUMBER PERFECT GAS FLOW 
3.1 Introduction 

Pop a hypersonic flow over a slender body at lower altitudes where 
density, and hence, Reynolds numbers are high, the flow will become 
turbulent. The Reynolds stresses and turbulent heat flux should be 
considered in the analysis of such flowfields. These two effects 
dominate the surface properties. However, at present, it is impossible 
to relate these fluctuating terms correctly to the dependent variables 
in the equations. Direct numerical solutions of turbulent flows cannot 
be obtained without a proper modeling of the ‘fluctuating terms. 

Many turbulence models have been developed with varying degrees of 
complexity [33]. These models, generally, are developed by first 
postulating a mathematical model containing undetermined constants, and 
then by attempting to choose the constants to make predictions fit the 
experimental measurements. Empirical turbulence models, such as the 
algebraic eddy viscosity models are appealing because the storage 
capacity and computing time required are much less than that for a more 
sophisticated turbulence model. Also, these models provide results 
which are comparable to a more complex model [33]. 

Two algebraic turbulence models, Cebeci-Smith [18] and Baldwin- 
Lomax [20], have been used widely for calculation of the compressible 
turbulent flows [19 ,22, 25, 3^-36] . Both are two-layer eddy-viscosity 



models and have similar forms. The primary difference between them is 
the choice of length and velocity scales in the outer layer. 

The boundary layer displacement thickness and boundary layer edge 
velocity are used as length and velocity scales, respectively, in the 
Cebeci-Smith model. The determination of the boundary layer edge is 
required within the solution. However, it is difficult to define the 
boundary layer edge in a hypersonic flow because there may not be a 
constant velocity region in the shock layer. Thus, the definition of 
the boundary layer edge is not well defined. Anderson and Moss used the 
ratio of the local total enthalpy to freestream total enthalpy to define 
the boundary layer edge [19]. It is based on the fact that the total 
enthalpy is constant within the inviscid part of the shock layer. It 
has been observed that with the presence of small numerical variations 
of the velocity or enthalpy profile in the shock layer, there could be 
sudden jumps and oscillations in the boundary layer thickness [37]. The 
length scale in the Cebeci-Smith model would be affected, hence the heat 
transfer to the body surface would experience oscillations that were not 
physically correct. 

Another definition suggested by Anderson and Moss [25] is based on 
the ratio of the integral of the viscous dissipation term in the shock 
layer for situations where the total enthalpy is not constant in the 
inviscid region. However, the factor of the ratio is subjective. 

By observing the shape of the total enthalpy profiles, Thompson et 
al. [3 1 *] defined the boundary layer thickness on the same physical basis 
as of Anderson and Moss [19] but in terms of the total enthalpy 
gradient. It was found that this definition gave reliable boundary 
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layer thickness consistently. It was superior to other definitions of 
the boundary layer edge. 

It has been shown that the calculated heat transfer was sensitive 
to the definition used in establishing the boundary layer thickness 
[373. Due to these difficulties with the Cebeci— Smith model, it is more 
appealing to use the Baldwin-Lomax model to predict the turbulent 
effects, since this model does not require the determination of the 
boundary layer edge. The velocity and length scales are based on the 
distribution of the vorticity. The maximum value of a vorticity 
function and its corresponding location, instead of the boundary layer 
quantities, are used to form the velocity and length scales for the 
outer layer. This model has become a popular algebraic eddy viscosity 

model. 

Although the Baldwin-Lomax model avoids determining the boundary 
layer edge in the outer eddy viscosity formulation, several difficulties 
have been encountered. First of all, there is ambiguity in determining 
the peak of the vorticity function. Degani and Schiff [38] found that 
there might be more than one peak in the function. This can result in 
an incorrect determination of the length scale for the outer layer if a 
wrong maximum value has been picked up. This problem can be eliminated 
by selecting the peak near the body surface [22,36]. 

The second difficulty lies in determining the two additional 
constants, C and C... . in the outer formulation. Baldwin and Lomax 

f Cp 

determined these constants by comparing results with the Cebeci-Smith 
model for transonic, constant pressure boundary layer flows [20]. 

Visbal and Knight [ 36 ] have shown that C cp should be decreased and C Kleb 
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increased for an equilibrium incompressible turbulent boundary layer 

[ 36 ]. However, Knight suggested a higher value of C for a compression 

cp 

corner calculation at Mach 3.0 [21]. The values of C and C V1 . depend 

cp Kl eb 

on the flow Mach number [20,21,36,37]. There is no one fixed value for 
all flow conditions. These values should be chosen carefully, otherwise 
different heat transfer results will be predicted [37]. These 
difficulties need to be investigated intensively before this model can 
be relied on to hypersonic flow conditions. 

In this chapter, a method for solving the flow over a blunt 
slender body where the inviscid region encompasses a significant portion 
of the total shock layer thickness is presented. The first order 
continuity and normal momentum equations are solved simultaneously as a 
coupled set rather than in a successive manner as has been utilized for 
wide-angle bodies. Two of the most frequently employed algebraic 
turbulence models, namely the Cebeci-Smith and Baldwin-Lomax models, are 
implemented to examine their suitability for the hypersonic flow. 

3.2 Basic Formulation 

As indicated in Chap. 2, the steady perfect gas viscous shock- 
layer equations [12] for an axisymmetric or two-dimensional body at zero 
angle of attack are obtained from the compressible Navier-Stokes 
equations, written in terms of a body-oriented coordinate system. They 
are nondimensionalized by the variables which are of order one in the 
region near the body surface (boundary layer) for large Reynolds 
numbers. The same set of equations are then nondimensionalized by 
variables which are of order one in the essentially inviscid region 
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outside the boundary layer. Terms in each set of equations up to 
second-order in the inverse square root of the Reynolds number are kept. 
A comparison of the two sets of equations is then made and one set of 
equations is obtained from them which is valid to second order in both 
the outer and inner regions. A solution of this set of equations is 
thus uniformly valid to second-order in the entire shock layer for 
arbitrary "Y. Anderson and Moss [193 used the eddy viscosity 
approximation to replace the Reynolds stresses and turbulent heat flux 
to find the solution of the turbulent flow. These equations provided 
here again are in an orthogonal, body oriented, transformed coordinates 
form, Eqs. (2.25) through (2.37). 

The second-order partial differential equations applicable to this 
study are expressed as [19] 

afw d 2 g/dn 2 * a, (dg/dn) « 2 

3n 2 (dg/dn) 2 3n (dg/dn) 2 


a 3 “4 3W 

+ si - — + - — -i- 

(dg/dn) 2 (dg/dn) 2 ^ 


(3-D 


where dg/dn and d 2 g/dn 2 are the first and second derivatives of the 
stretching function g. The quantity W represents u for the s-momentum 
equation and H for the energy equation. The coefficients to a 1) are 
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written as: 


s-momentum, W - u: 


1 

y(1+e + ) 


3y(1+e ) dg + 
3n dff 


d + n n sh K) (1 + e + ) 


jn .cose 
sh 

r + Tin .cose 
sn 


n' n . npu 
sh sh r 


2 + — 
e y(1+e )(1+nn ic) 


n . pv 
sir 

2 + 
e y(1+e ) 


(3-2) 


"sh* 


p(1+e )(1+nn gh K) 


3y(1+e ) dg 
3n dff 


2 

n . k 
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(1+e )(1+nn gh K) 


jcose 


1+nh sh K 


r+fih .cose 
sh 


) - 


n sh Kpv 


2 — + 

e ( 1+ n n ah <)y( 1+ e ) 


(3.3) 



2 — + 
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sir 


H 2 — + 
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Energy, W « H: 
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The remaining first-order equations are written as 
Global continuity: 

f^ Cn 3h (r+T,n sh c03e)Jpu:i + tfT 3^ {(r+nn sh co30) 


( 3 . 6 ) 


( 3 . 7 ) 


( 3 . 8 ) 


( 3 . 9 ) 


( 3 - 10 ) 
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x [( 1 +nn ah <)pv - n^ h npu]} - 0 


(3.11) 


n-momentum: 


pu 9v 
(Unn^) 9£ 


n tfi npu 9v dg + 
n sh (1+ ^sh K) 3ti d7T 


pv 9v dg 
n 3h 3n ^ 


2 

PU K 


(1+rih 



+ J_ 9£ dg 
n 3h 3ti dTf 


0 


( 3 . 12 ) 


State: 

p - pT Or-n/Y (3.13) 

The molecular viscosity is given by the Sutherland's law as 

U - [(1+C)/(T+C)]T 3/2 (3.14) 

where 

C-cVdhV (3.15) 

and C is 110.33 K for air. 

In the preceding equations, the prime denotes the differentiation 

with respect to £, and e + is the eddy viscosity which is set equal to 
zero for a laminar flow. The independent variable of transformation is 
defined by 

n - g(rp-) - g(n) (3.16) 

sh 

The stretching function g(n) , Eq. (2.21), is given by (a - 0) 



36 


g(n) “ 1 “ 


1 .-r B-h+1 

--ft: 

ln( Fi > 


and it3 first and second derivatives are 

1 r 1 j. 1 1 

t rm i ™. rr + 7TFr®nrrT J 


d£ - 
dn 


TTFiT 1 (I- n+i~) + (? + n-i )" 


(3.17) 


( 3 . 18 ) 


d^g . 1 [ 1 _ - _J ] ( 3 . 19 ) 

dn 2 (6-n+1) 2 (e + n~1) 

Equation (3.17) permits the mesh to be refined near the body with 
the values of I near 1 giving the largest amount of stretching. 

Equation (3.16) may be inverted to obtain the physical coordinate n from 
the transformed coordinate n: 


- r l FT J 
n - 1 - B [ — Z 


— „ 1-n 
eili _ i 


„ . 1 ~n 

<£f> *’ 


( 3 - 20 ) 


The transformation of Eq. (3-16) keeps the body at n - 0 and the shock 
at n - 1 with uniform mesh in the computational coordinate n. 


3.3 Boundary Conditions 

At low altitude, slip effects are not important. At the wall, no' 
slip and no-temperature-jump boundary conditions are used. The wall 
temperature and enthalpy are specified as constant. 
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The boundary conditions at the shock are calculated by using the 
Rankine-Hugoniot relations. The nondimensional forms are given as 
Continuity: 

"sh ? 3h ■ ' 3l " a (3 - 2,) 

T angent i al -moment um : 

Q - cosa (3-22) 

sh 

N ormal -momentum: 


P 


sh 


— — r + sin 2 a (1 
YM 2 

ao 



(3.23) 


State 


sh 


^ Y ^ p sh T sh 


Density 


sh 


(Y+1 )M 2 sin 2 a 

GO 

(Y-1 )M 2 sin 2 a + 2 


(3.24) 


(3.25) 


where Q g ^ and are the velocity components in the tangential and 

normal directions, respectively, in the shock-oriented coordination 
system. These are related to the body-oriented coordinate as 

u sh “ Q sh sin (“ + 0) + v 3h cos(a+B) (3.26) 


- a gh cos(a+B) + V gh sin(a+B) 


(3.27) 
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3.4 Turbulence Models 

Two of the most widely used algebraic turbulence models, namely 
the Cebeci-Smith (C S) and Baldwin-Lomax (BL), have been implemented in 
this study. Algebraic turbulence models are more appealing because they 
require less computer storage and much less computational time as 
compared to the two-equation model of turbulence, such as <-e model. 

Both Cebeci-Smith and Baldwin-Lomax models of turbulence employ a two- 
layer eddy-viscosity formulation. The inner law is based upon Prandtl's 
mixing-length concept. The outer law employs either the Clauser- 
Klebanoff expression (in the Cebeci-Smith model) or an equivalent 
expression (in the Baldwin-Lomax model) for computing the eddy 
viscosity. 

3 . 4 .I Cebeci-Smith Turbulence Model 

The algebraic eddy viscosity (in nondimensional form) is given by 

+ — ^ — 

e i n n crossover 

( 3 . 28 ) 

+ — — 

£ o n ^crossover 

where n is the value of n at which values from the inner and 

'crossover 

outer formulas are equal. 

The inner eddy viscosity is obtained from the Prandtl mixing- 
length concept 




( 3 . 29 ) 
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The mixing length & is obtained by using the van Driest' s proposal 
stated as [18,19] 

A - K i n 3h n[1 - exp(-n + /A + )] (3-30) 


where 


n 


n sh np 

eu 


w 


3u dg 

P n sh ^ ^ ^ j w ] 


1/2 


(3.31 ) 


The quantity is the von Karman constant with a value of 0.M, and A 

is a damping factor expressed (for flows with a pressure gradient) as 
[18,19] 


A + - 26(1-1 1 .8P + ) 1/2 


where 


t e sh e e 


and 


y 3u dg 1/2 

„ . e [ — *£_ (_ _) ] 

T Pn 3h 3,1 « 


(3.32) 


(3.33) 


(3.3*0 


For the outer region of the viscous layer, the eddy viscosity is 
approximated by the Clauser-Klebanoff expression [18,19] 

+ 
e 

o 

where 


K pu 5. Y. - 
2 e k i,n 
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e U 
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(3.36) 


0.0168 


(3.37) 
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and 


i » n 


n . n 

►5.5(-Sji-) 1 


6 -1 


(3.38) 


The boundary layer thickness, 6, in Eqs, (3*36) and (3*38) is 

assumed to be the value of ti at the point where 

H. / H. - 0.995 (3 ' 39 

t t ,°° 

Another criterion for obtaining 6 , based on viscous dissipation, is the 
height where 


r° t 

J 0 (dg/dn) 

Z 1 . X — 

J 0 (dg/dn) 


dn 


dn 


0.995 


(3.^0) 


where 


F - { e 2 y [(1 + e + )~- fi* - 
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icu 
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sh 
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2 /(1+e + )u 


(3.^1 ) 


Thompson et al. [3*0 defined the boundary layer edge at a location where 
d (H t / H t J/dn S 0.5 (3 ‘ 42) 

All three of these criteria have been used in this study. 

3.4.2 Baldwin-Lomax Turbulence Model 

This model employs a formulation similar to the Cebeci-Smith model 

for the inner-region eddy viscosity 

«: - M (3 - 1,3) 

e u 

where l is given by Eq. (3-30) except that A + is replaced by A + which is 


defined as 
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In Eq. (3. MM), t is the local shear stress obtained from 
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The magnitude of the vorticity, |u>| is given by 


1 3v qn^ h dg 3v 

^ I (1+nn sh K)^95 n„ u dn 9n 
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(1+n n gh<) dg 3u 

= <u] 

sh dn 3n 


The outer-eddy-viscosity approximation of the Baldwin-Lomax 
replaces the Clauser-Klebanoff formulation by the relation 

+ K 2 C cp pF wake F KLEB (n) 

£ 0 2 
e y 


where is a constant given by Eq. (3.37), C cp is an additional 

constant given as 1 .6 [20] and 

F ■ n F 
wake 'max max 

The quantities n m _„ and F are the values at the location of the 

UlaX ulaX 

maximum value in the vorticity function 

F(n) - n 3h n| <»)| (l-exp(-n + /A + ) ) 

The Klebanoff intermittency factor, F , is given by 

KLEB 


F KLEB (n) 


(1+5. 5( 
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" here C Kleb 


0.3 


It has been found that the value of C depends on the Mach 

cp 

number. Baldwin and Lomax [20] chose a value of 1.6 for C cp by 

comparing with the Cebeci-Smith' s turbulence model for transonic, 
constant pressure boundary layer flows. However, a value of 3.0 for C cp 


is found more appropriate for hypersonic flow [37]. 

3 . 4.3 Transition Model 

Both continuous and instantaneous transitions from a laminar to 
turbulent flow have been included in this study. Instantaneous 
transition is initialized when the local Reynolds number or momentum- 
thickness Reynolds number exceeds a preselected value. Continuous 
transition is effected by defining a streamwise transition intermittency 

+ 

factor Y, which modifies the composite eddy viscosity e over a 
1 » 5 

specific distance along the body. 

The factor Y, is evaluated by a relation developed by Dhawan and 
* • £ 

Narasimha [39] as 

Y - 1-exp(-0.4l2O (3.51 ) 


where 


- *<e-V 
5 * e 0 «-u 


(3.52) 


The quantity is the location where the transition is started and x 


is approximately equal to two. 
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3.5 Method of Solution 

The overall method of solution employed is an implicit finite- 
difference spatial-marching method, similar to the one employed by 
Davis, and Anderson and Moss [12,19]. However, the method is 
implemented in this study differently because convergence problems are 
encountered for slender bodies if the method of solution outlined in 
these references is employed. 

To simplify the numerical computations, the viscous shock-layer 
equations are again transformed by normalizing most of the variables 
with their local shock values. It should be pointed out that the normal 
velocity, at the shock may change sign at some locations and may be near 
zero at others. The normalized v-profiles in such a region are not very 
well behaved and stability problems can occur if these profiles are used 
in the solution procedure. Therefore, it is desirable to remove the 
normalization procedure from the normal velocity profile. 

When the normal coordinate is normalized with respect to the local 
shock standoff distance, a constant number of finite-difference grid 
points between the body and shock are used. The second order equations, 
Eqs. (3.1) through (3*10), are solved by using the finite difference 
method. The derivatives are replaced with finite-difference expressions 
in a such way that three-point central differences in the n~di^ection 
and two-point backward differences in the ^-direction occur. The 
truncation terms of order A5 m and either ^n n An n _ 1 or A h n _ An n _ 1 are 

neglected. The subscript n denotes the grid point along a line normal 
to the body surface, whereas the subscript m denotes the grid station 
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along the body surface. Replacing the differential terms by the finite - 
difference expressions, the governing equations are expressed as 


A „Vn-1 * V.,n * C nVn*1 ’ D n 


(3.53) 


where 


^ a 1 ^m,n Ari n 


n " 4 Vt ( ‘V‘Vi > 4 Vi ( ‘V*Vi 1 


(3.54) 
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(3.56) 


D 

n 


(3.57) 


Equation (3.53) along with the boundary conditions constitutes a system 
of the tridiagonal form, for which efficient computational procedures 
are available. 

To avoid the instability problem encountered by the traditional 
approach [12,19] and from the hypersonic small disturbance theory [14] 
for a flow on a slender body, the two first-order equations — the 
continuity and normal-momentum equations — are solved simultaneously as 
a coupled set rather than in a successive manner for the pressure, p, 
and normal velocity, v. The density in these equations is eliminated 
through use of the equation of state. The resulting equations are 
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(3.59) 


Equations (3.58) and (3.59) are expressed in the finite-difference form 
at points (m,n+1/2) and (m,n-1/2) using a box scheme discussed by 
Richtmyer [MO]. The final forms are 


A c,n+1/2 v m,n+1 

+ D , p 
c,n+1/2 F m,n 


B c, n+1/2 v m,n + C c, n+1/2 p m,n+1 
E c f n+1/2 


(3.60) 
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( 3 . 61 ) 
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(3-63) 


The coefficients of these equations are given in Appendix A. Eliminating 
p and v alternatively in the coupled equations, Eqs. (3.60)-(3.63) « two 
tridiagonal equations for pressure and normal velocity are obtained as 


m,n+1 

+ K 2 

P + 

K m,n 

K 3 

p m, 

n-1 


+ K, 

v + 

K„ 

V 


m,n+1 

6 

m,n 

7 

m, 

n-1 


The coefficients K 1 through K g are also given in Appendix A. 


(3.64) 

(3-65) 

Equations 


(3 # 64) and (3.65) are solved in the same way as the energy and s- 


momentum equations. 

By integrating the continuity equation from n-0 to n-1 , a 
quadratic equation for the shock standoff distance is obtained. The 
density is determined by the equation of state. 

The solution is started at the stagnation streamline where the 
various flowfield quantities are expanded in terms of the distance, 5, 
along the body surface [19,26]. These series expansions reduce the 
partial differential equations, Eqs. (3.1) through (3.12), to ordinary 
differential equations in terms of n. At a body location m, other than 
the stagnation streamline, a two-point backward difference scheme is 
used for the derivative with respect to £ at the point (m,n). This 
again gives ordinary differential equations at location m in terms of n 
for Eqs. (3.1) through (3.12). The finite-difference form of these 
ordinary differential equations (obtained through the central 
differences) can be solved by using the Thomas Algorithm. Figure 3.1 
gives the flow chart for the solution sequence of these equations. 
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The Vigneron condition [41] (for the pressure gradient in the 
streamwise momentum equation) has been used for marching in the subsonic 
nose region. In this condition, a portion of the pressure gradient is 
treated implicitly by employing a two-point backward difference scheme. 
The remainder portion of the pressure gradient is forward differenced to 
allow for upstream influences. 

The solution is iterated at location m until convergence is 
achieved. The solution advances to the next body station, m+1, and uses 
the previous converged solution profiles as initial values for starting 
the solution at m+1 . This procedure is repeated until a global solution 
at all body locations is obtained. 

The initial shock shape is created by the thin layer approximation 

0 

for a short wide-angle body (35 sphere-cone, for example). The shock 
shape obtained from a full-layer solution to this body shape is then 
used as an initial guess for the slender bodies in a sequential manner 

by reducing the body angles in steps of 5 to 10 degrees. In place of 
the shock stand-off distance used previously [17,42], its derivative in 
the streamwise direction is smoothed after each global iteration. 

Due to the change in sign of the normal velocity profile from 
station to station, an under-relaxation scheme [43] 

F - 55F 1 + (1-w) F 2 (3-66) 

has been employed in the present work. Here F 1 is the most recently 

calculated physical quantity and F 2 is the value obtained from the 

previous local iteration. A value of 2 of 0.2 to 0.4 gives convergence 
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in most cases. In general, an under-relaxation was required only for 
the pressure and normal velocity. 

Depending upon the initial approximation to the shock shape 
(whether obtained by using the thin shock-layer form of normal -momentum 
equation or from a larger body angle solution) , the first global pass 
solution may be improved by subsequent global iterations. 

3.6 Results and Discussions 

Numerical solutions to the previously discussed viscous shock- 
layer equations for the hypersonic flow over a long slender body have 
been obtained. Results for laminar, transitional, and turbulent flows 
of a perfect gas are compared with the experimental data and/or with 

numerical solutions in the literature. The solutions are chosen for 

0 0 

small body angle (5 to 35 ) hyperboloids and spher.e-cones at zero- 
degree angle of attack. The free stream Reynolds numbers are within the 

range from 1.2 X 10 4 to 3-5 X 10 6 . 

3.6.1 Comparison of the Present Method with Cascading Method 

Figures 3.2 to 3.6 show the effect of solving the normal momentum 
and continuity equations simultaneously in a coupled way as compared to 
solving all the governing equations in a successive way [12,19]. 

Results of shock stand-off distance, wall pressure, and skin friction 

0 

coefficient are shown in Figs. 3.2 to 3. 1 * for a hyperboloid with 20 
half-body angle. It is clearly noticed that the solutions oscillate in 
the downstream region with the cascading approach. It is also noticed 
that this instability can be removed when coupling is implemented 
between the normal momentum and continuity equations. 
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Figures 3.5 and 3.6 show the results for shock stand-off distance 

0 

and Stanton number distribution for a 35 sphere-cone without coupling 
the two first-order equations. Oscillation in the solution is noted in 
the vicinity of the tangency point where the curvature is discontinuous. 
The curvature is equal to one on the spherical part and zero on the 
conical part. It is important to note that coupling the two first-order 
equations can stabilize the solutions at this discontinuity. 

3.6.2 Comparison of the Present Method with other Predictions and Data 
for Laminar flows 

The results obtained by the present method (VSL2D) are compared 
with another method (VSL3D) in Figs. 3.7 to 3.10. The results of VSL3D 
were obtained by Thompson [44]. Figure 3*7 gives the convergence 
history of the streamwise derivative of the shock stand-off distance for 
the present method. It is seen that solutions do not diverge with the 
subsequent global passes. Figure 3.8 presents a comparison of the 
boundary layer thickness, as obtained by the VSL2D and VSL3D methods. 

The boundary layer thickness is defined as the location where “ 

0.995. The predictions of the boundary layer thickness are quite 
different by the two methods. The VSL3D results show a big jump in the 
boundary layer thickness. In spite of this difference, the laminar 
heat transfer and skin-friction coefficients compare well as shown in 
Figs. 3-9 and 3-10. 

0 

A comparison of surface heat transfer results for a long 5 
sphere-cone between the present viscous shock-layer (VSL) and 
parabolized Navier-Stokes (PNS) predictions [45] is given in Fig. 3.11. 
The VSL results are about five to ten percent higher for most of the 



20* hyperboloid 
With coupling 
Without coupling 
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Flo. 3.3 Wall pressure distribution for a 20* hyperboloid. 



20* hyperboloid 
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Flo. 3.4 Skin-friction coefficient for a 20* hyperboloid. 
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Fig. 36 Stanton number distribution for 
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Fig. 3.9 Surface heat transfer rate distribution for 
sphere-cone with laminar flow. 
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Fig. 3.10 Skin-friction coefficient for a io* sphere-cone 
with laminar flow. 
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body length as compared to the PNS results. The PNS predictions employ 
fourth-order explicit and second-order implicit smoothing terms, whereas 
the present VSL calculations do not use any smoothing. Furthermore, the 
stability of PNS solutions restricts the reduction of the normal grid 
spacing adjacent to the wall (required for accurate heat transfer 
predictions) if a relatively large marching stepsize is required for a 
long body. Since the PNS requires a starting solution that describes 
the subsonic region, any starting solution errors distort the PNS 
results in the nose region. In the VSL method, the starting profiles 
are created as part of the solution and, thus, the method is self- 
starting. 

Figures 3.12 to 3.15 show comparisons of results obtained from the 

present method with available experimental data. A comparison of the 

0 

predicted pressure distributions on a 10 hyperboloid -with the 
experimental data [46], as well as with the results of Hosny et al. 

[15], is given in Fig. 3.12. Both predictions compare quite well with 
the data. Hosny et al. [15] solved all governing equations in a coupled 
manner which may require more computational time at every point in the 
flow, especially if real gas properties are included. The present 
method with the coupling between the two first-order equations gives 
equally accurate results. The present approach may be more appealing 
for real gas calculations where local iterations are required to update 
the chemical composition along with the transport and thermodynamic 
properties. 

The present method gives surface heat transfer rates which compare 

0 

quite well with the data of Cleary [47] for a 15 sphere-cone as shown 
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with laminar flow. 




Fig. 3.14 Surface pressure distribution for a 12.84 
sphere-cone with laminar flow. 
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in Fig. 3.13. As shown in Figs. 3.14 and 3-15, the present predictions 
for the wall pressure distribution and surface heat transfer rate on a 

12.84° sphere-cone agree fairly well with the experimental values of 
Miller [48,49]. 

3 . 6.3 Comparison of Predictions with different Turbulence Models and 
Data for Long Slender Bodies 

Cebeci-Smith [18] and Baldwin-Lomax [20] turbulence models are 
implemented in the present method to predict turbulence effects. 
Transition to turbulence is modeled by using the Dhawan and Narasimha 

method [39]. 

Figures 3 . 16 to 3.18 show a comparison between the results for a 

10 sphere-cone as obtained by the VSL2D and VSL3D [44] models. The 
Cebeci-Smith turbulence model is implemented in both solutions. The 

onset of transition is set at s /R^ • 2.0. The definition of the 

boundary layer edge is based on the total enthalpy (i.e., * 

0.995). The results for the boundary layer thickness are shown in Fig. 
3.16. Similar to the results for the laminar flow (Fig. 3.8), these 
results also differ significantly with each other. However, this 
difference influences the surface heat transfer rate and skin friction 
coefficient predictions significantly for a turbulent flow as shown in 
Figs. 3.17 and 3-18. 

The length and velocity scales in the Cebeci-Smith model are 
strongly dependent on the boundary layer edge location, and these 
influence the surface properties. Due to this difficulty, Thompson et 
al. [34] defined the boundary-layer edge location based on the gradient 
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of the total enthalpy [ 3 (H /H ) /3n] . With this definition, good 

v U | 

agreement between the VSL2D and VSL3D results is obtained for the 
surface heat transfer rate (Fig. 3.19). The results obtained by using 
the classical total enthalpy definition for the boundary-layer edge are 
also shown in Fig. 3.19. It is noted that the new boundary layer 
definition gives results comparable to the classical definition along a 
long body [37]. 

0 

The results for the turbulent flow over a 9 sphere-cone are 
illustrated in Figs. 3.20 to 3.25. The hemispherical portion of the 
model was roughened in order to insure attainment of turbulence flow 
over this region [50]. Three different definitions, which are based on 
the total enthalpy, the gradient of total enthalpy, and the dissipation 
models, for the boundary layer edge locations have been used to 
calculate the turbulent heating with the Cebeci-Smith model as given in 
Fig. 3.20. It is seen clearly that predictions from the total enthalpy 
and total enthalpy gradient models are comparable to the experimental 
data. The boundary layer thickness based on the dissipation model does 
not give heat transfer predictions comparable to the classical enthalpy 
model . 

The Baldwin-Lomax model uses the distribution of vorticity to form 
the outer length scale. It is known that there may be more than one 
peak in the vorticity function [38]. Figures 3-21 to 3.23 give the 
distributions of the vorticity function at three different locations. 

It is seen that there is more than one peak at each location. The 
correct peak that should be picked is near the body surface. The 
Baldwin-Lomax model proposed originally [20] was for constant pressure 
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boundary layers at transonic speeds. This model has been modified to 

include the effect of pressure gradient on the damping factor A [37]. 
Results of surface heat transfer for both models are given in Fig. 3-24. 
When compared to the experimental data, the original model predicts well 

up to the tangency point (s*/R N *« 1.5) in the favorable pressure 

gradient region, whereas the modified model gives good predictions in 

the adverse pr essure _ gradient region and beyond. It was, therefore, 

decided to combine the two models. The predictions for the combined 

pressure gradient models implemented in the Baldwin— Lomax model are 

given in Fig. 3.25. Here the original Baldwin-Lomax model is used up to 

the tangency point and the modified model is used afterwards. The 

combined model gives very good predictions when compared to the 

experimental data along the entire body length. Also included in this 

figure are the predictions obtained with the Cebeci— Smith model. The 

two models give almost the same surface heat transfer predictions. 

0 

The results for the Stanton number distributions for a 7 sphere- 
cone are given in Figs. 3.26 and 3-27. The transition to turbulence is 

initialized at s /R. t - 4.8 as given by the data of Carver [51]. 

Results by using the Cebeci— Smith turbulence model with two boundary 
layer edge definitions based on the total enthalpy and its gradient are 
shown in Fig. 3.26. Both definitions give surface heat transfer rate 
predictions within 15$ to the experimental data. There is an increase 

ft * 

in the value of Stanton number at a s /Rjj location of about 15 by using 

the total enthalpy definition for the boundary layer thickness. This is 
probably due to the poor resolution of the gradients of various 
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Fig. 3.22 Distribution of vortlclty function at s - 2.0 for 
a 9* sphere-cone with turbulent flow. 
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Fig. 



F<n> 

3 23 Distribution of vorticity function at s • 10. 0 for 
a 9» sphere-cone with turbulent flow. 
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flowfield quantities at the boundary layer edge at this body location. 

This problem can be overcome through the use of adaptive grids. The 

observed increase in the heating data in Figs. 3*26 and 3.27 with 

increasing body location is probably due to the entropy-layer 

swallowing. This trend is also visible from the theoretical 

predictions. Once the swallowing is complete, the data and predictions 

would probably show a decreasing trend with increasing body distance. 

Since the sensitivity of the surface heat transfer rate on the 

boundary layer thickness definition does not exist in the Baldwin-Lomax 

turbulence model, the Stanton number distribution was obtained using 

this model in Fig. 3.27. A comparison of predictions with the 

experimental data shows that a value of 3.0 for C gives better results 

cp 

as compared to the values of 2.08 suggested for * 3 by Knight [21] or 

the original value of 1.6 given for M m - 1 [20]. This coefficient 

which appears in the Baldwin-Lomax outer formulation is dependent upon 

the flow Mach number. The value of 3 inferred here for C at M » 8 

cp “ 

along with the other suggested values for different Mach numbers point 
to a linear dependence of C cp on the flow Mach number in the range 1 £ M 

£ 8. Additional comparisons with data are necessary to verify this 
dependence. 

0 

The surface heat transfer results for a 5 sphere-cone obtained by 
using the Cebeci-Smith and the modified Baldwin-Lomax models of 
turbulence are illustrated in Fig. 3.28. The results of laminar flow 
calculations are shown in this figure. These calculations were 
performed using both the present method and the VSL3D method [3^]. The 

C-2 





Fig. 3.26 Stanton number distribution for a 7 * sphere-cone 
with turbulent flow. 
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Flo* 3.27 Stanton number distribution for a 7* sphere-cone 
with turbulent flow. 
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transition was initialized at an axial location of 23.8 ft. Even though 
the turbulent flow results from the two methods using the Cebeci-Smith 
model are different, the VSL3D results compare with the modified 
Baldwin-Lomax model very well from the end of transition. 

3.7 Conclusions 

Numerical solutions of viscous shock-layer equations are presented 
for hypersonic laminar and turbulent flows over long slender bodies. 
These results are obtained from a method which employs a spatial- 
marching implicit finite-difference technique. This technique is fast 
and uses partial coupling among the governing equations based on the 
hypersonic small disturbance theory. The partial coupling yields a 
simple and computationally efficient technique. 

Detailed comparisons have been made with other predictions and 
experimental data for slender body flows to assess the accuracy of the 
present numerical technique. Results from the present method show that 
the coupling between the normal momentum and continuity equations is 
essential and adequate to obtain stable and fairly accurate solutions 
past long slender bodies. 

The two widely used algebraic turbulence models, namely, the 
Cebeci-Smith and Baldwin-Lomax models have been analyzed with the 
present numerical technique for application to long slender bodies. 

Both of these models appear adequate for such flows. Due to the 
sensitivity of the Cebeci-Smith turbulence model to the boundary layer 
edge location, however, it is imperative that the numerical method 
should provide good resolution and accurate solutions near the boundary 
layer edge. This can be a problem for long slender bodies, especially, 
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if the numerical method (such as PNS) employs artificial viscosity to 
damp oscillations. For this reason, the Baldwin-Lomax turbulence model, 
which avoids the use of the conventional boundary layer thickness in its 
formulation, appears more convenient to implement. 

A correction for the pressure-gradient effect has been made to the 
Baldwin-Lomax model. Constant C in the outer-layer formulation has 

cp 

been modified to 3.0 for the Mach 8 case. Based upon this study and 
other investigations, a linear dependence of this constant on the flow 
Mach number is suggested. Further comparisons with the experimental 
data are needed to verify this dependence. An additional consideration 
in the implementation of the Baldwin-Lomax model concerns the appearance 
of two peaks associated with the two maxima in the vorticity functions 
used to form the outer-layer length scale. The second peak is avoided 
by choosing the first one in the region where the gradient of total 
enthalpy is less than or equal to 0.5, i.e., 3(H t /H t ^ oo )/3ri £ 0.5. 


85 


Chapter 4 

LOW REYNOLDS NUMBER PERFECT GAS FLOW 
4.1 Introduction 

Most future hypersonic vehicles will be operating in the upper 
atmosphere, where "low density effects" will play a major role in 
establishing the lift, drag, moments, and aerodynamic heating on a 
hypersonic vehicle. An accurate knowledge of hypersonic 
aerothermodynamics under low density conditions is required for an 
accurate prediction of the aerothermal environment for the new 
generation of hypersonic vehicles. 

The degree of rarefaction of a low density flow is usually 
expressed through the Knudsen number which is the ratio of the molecular 
mean free path in the gas to a characteristic dimension of the 
flowfield. The conventional continuum flow assumption is valid when 
this parameter is very small in comparison to unity. The opposite limit 
of very large Knudsen number corresponds to a free molecule flow in 
which intermolecular collisions may be neglected. The region between 
these limits is generally referred to as the transition flow regime 
[52]. 

At highly rarefied gas flow conditions, the conventional approach 
by continuum analysis is no longer valid. A more appropriate approach 
is by kinetic theory of gases which can correctly describe microscopic 
properties of molecules, such as the Boltzmann equation [30]. Although 
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there is only one dependent variable, the distribution function for the 
molecular states in the Boltzmann equation, the number of independent 
variables make this approach extremely difficult to obtain analytical or 
numerical solutions. An alternative is to model the gas flow at the 
molecular level. The Direct Simulation Monte Carlo (DSMC) method [531 
has been found to be most readily applicable to complex engineering 
problems. However, this method still requires large computational times 
and computer storage. 

At the slightly rarefied gas conditions, significant levels of 
molecular collisions are still present in the flowfield which make the 
continuum approaches applicable except in a region next to the wall. It 
is because the gradients of the macroscopic variables become so steep 
that the mean free path becomes large compared to the local 
characteristic length. This region is called the Knudsen layer in which 
the determination of the flow properties requires the direct solution of 
the Boltzmann equation matched to the solutions for the outer flow and 
the wall boundary conditions. This is most conveniently done through 
the use of a slip model in which slip and jump properties are used for 
the boundary conditions for the conventional continuum flow equations. 
These slip and jump boundary conditions for the gas and solid interface 
are obtained from the balance equations for mass, momentum and energy 
fluxes at the Knudsen layer edge C 30—31 , 5M-55]. 

Not much attention seems to have been given to the problems 
encountered with low-density aerothermodynamics. Davis [12] included 
body and shock-slip in the viscous shock-layer analysis of a perfect-gas 
flow around a hyperboloid. Davis [56] modified these slip relations for 
a binary mixture. Tree et al. [57] analyzed the hypersonic ionizing 
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viscous shock-layer flow past axially symmetric bodies at low densities. 
Tiwari and Szema [58,59] investigated the effects of body and shock slip 
conditions on the aerothermodynamic environment of a Jovian entry body. 
Swaminathan et al. [60] and Song et al. [61] recently included the body- 
and shock-slip effects for three-dimensional flows. However, their 
surface slip condition for single species or multicomponent mixtures 
contained some errors as explained in Ref. 31- The shock-slip boundary 
conditions did not account for the derivatives of the shock quantities 
in the shock-oriented coordinate system. This introduced significant 
errors in analyzing flows past slender bodies as compared to the wide- 
angle bodies [62]. 

Gupta et al. [31] reanalysed the wall boundary conditions by using 
the approach of Scott [63] and provided appropriate relations for the 
various quantities with surface slip in a form which can readily be 
employed for multicomponent and binary mixtures as well as a single- 
species gas. These surface slip expressions have been implemented 
successfully in full Navier-Stokes equations [26,62] 

Under the low Reynolds number (or low density) flow conditions, 
the viscous effects influence almost the entire 3 hock layer and the 
shock itself is considerably thick as compared to the high Reynolds 
number (or high density) case. The complete Navier-Stokes equations are 
considered appropriate for the low Reynolds number applications. But 
computer storage and computational time make them are very expensive to 
solve for flows around long bodies. The viscous shock-layer equations 
have been shown to give good results for hypersonic flows on blunt 
slender long bodies at high Reynolds number. Thus, it is desirable to 
employ the viscous shock-layer equations instead of Navier-Stokes 
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equations at low Reynolds number cases to reduce the computational 
requirement. 

In this chapter, the surface-slip relations developed by Gupta et 
al. [31] and the corrected form of the shock-slip boundary conditions 
[12] are implemented in the viscous shock-layer code for a perfect gas 
as described in Chap. 2. to obtain results for the low-density flight 
conditions for long slender bodies. A detailed comparison with 
experimental data and other numerical results gives an estimate of 
accuracy of the present predictions. Furthermore, the range of 
applicability of the viscous shock-layer solutions is ascertained by 
comparing these results with those obtained from the steady-state 
Navier-Stokes equations. 

4.2 Flow Governing Equations 

The conservation equations employed in this chapter are the steady 
perfect gas viscous shock-layer equations for an axisymmetric or two- 
dimensional body at zero angle of attack [12]. These equations are 
written in the same forms as Eqs. (2.25) through (2.37) except that the 

eddy viscosity, e + , is set to zero. Also, the stretching function g(n), 
Eq. (2.21), is given by 


(1-q) 

an( FT ) 


An { 


B~n( 2a+1 )+1 

F+nX 2a+1 )-1 


}] 


g(n) “ 1 - [a 


(4.1) 
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The first and second derivatives of Eq. (4.1) are expressed as 


dg (1-ot) (2a+1 ) r 1 i 

drf * „_, lf+1 ^ 1 [^— n"( 2a+1 ) + 1 ] [£+n(2a+1 )-1 ] i 

z ( I r T ) 


d 2 g _ (1-q)(2a+1) 2 | 1 1 

dn 2 Jtn(|^j) [e-n(2a+1)+1] 2 [B+n(2a+1 )-l ] 2 


Equation (4.1) permits the mesh to be refined either near the body 
only (a - 0) or refined equally near both the body and bow shock (a - 
1/2) when the shock becomes thick under the low density flight 
conditions. The parameter B controls the amount of refinement with 
values near 1 giving the largest amount of stretching. The physical 
coordinate n can be obtained from inverting Eq. (4.1) and is expressed 
as 


(25+1) 
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(4.4) 


+ 1 


This transformation keeps the body at n - 0 and the shock at n ■ 1 with 
uniform mesh in the computational coordinate n. 
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4.3 Shock- and Surface-Slip Boundary Conditions 
At high altitudes, the continuum flow assumption breaks down in 
the region next to the wall. The no-slip and no-temperature jump 
boundary conditions are no longer valid. As such, the slip and 
temperature jump boundary conditions should be used. The relations for 
the body and shock slip conditions are provided in this section. 

4.3.1 Surface Slip Conditions 

The surface slip conditions for a single-species gas as given in 
Gupta et al. [31] are used as the boundary conditions on the body 
surface. Since no mass injection is considered in this chapter, the 
normal component of velocity at the surface is taken to be zero. The 
nond iraensional forms of surface slip conditions in the computational 
plane, Eqs. (2.59) through (2.61), are given here again as 
Velocity slip: 


u 

3 



[ 


1 dg 3u 
n sh dp 3n ’ 


KU 




(4.5) 


Pressure slip: 
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Temperature slip: 
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In the derivation of the relation for the temps' ature slip, it is 
assumed that the internal energy is frozen during the reflection from 
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the wall. The parameter 0 is the accommodation coefficient which is 
taken to be 1 in this study. 

4. 3.2 Shock Slip Conditions 

The boundary conditions at the shock are the modified Rankine- 
Hugoniot relations developed by Cheng [32]. These relations are 
obtained by integrating the one-dimensional Navier-Stokes equation 
across the shock. Since velocity components tangent and normal to the 
shock are not the same as those tangent and normal to the body surface, 
transformations are needed to relate these quantites. The 
transformations are given by 

Q sh “ u sh CO3 ( a-0 ) + v gh sin(o-0) (4.8) 

7 sh “ " u sh sln(a-9) + v gh cos(a-0) (4.9) 

where Q gh and ? gh are the components of velocity tangent and normal to 

the shock interface, respectively. The nondimensional forms of the 
shock slip conditions in the computational plane, Eqs. (2.75) through 
(2.78), are given here again by 
Continuity: 


p . V . « -sina 
r sh sh 


Tangential-momentum: 


(4.10) 
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(4.11) 
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Normal-momentum: 
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Equation of state: 

»»h- y »3h /(Y -’ )T 3h 

The errors in the expressions U3ed by Davis [12] have been 
corrected in this study and this is discussed also by Lee et al. [62]. 

4.4 Method of Solution 

The method of solution is similar to that implemented in Chap. 3. 
The two second-order equations, s-momentum and energy, are replaced with 
central differences in the q-direction and two-point backward 
differences in the ^-direction. The two first-order equations, 
continuity and normal-momentum, are solved simultaneously in a coupled 
way. The solution is started at the stagnation point. The velocity 
slip and temperature jump on the surface and at the shock are iterated 
along with the corresponding governing equations. The solution is 
iterated at location m until the convergence is achieved. The solution 
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is then advanced to the m+1 station. Figure 4.1 presents the flow chart 
for obtaining solutions with body-slip only, shock-slip only, and with 
body and shock-slip. 

Solutions to the steady-state Navier-Stokes equations, which are 
given in Appendix B, have been obtained by first expressing them in the 
body-oriented coordinate system. This procedure is the same as the one 
employed with the viscous shock-layer equations. After obtaining a 
solution of the viscous shock-layer equation, the higher-order terms are 
evaluated using these flowfield results. These terms are held constant 
during the solution for the first approximation to the Navier-Stokes 
equations. The solution with this approximation is obtained at the end 
of the first global pass. At the beginning of the second global pass, 
the higher-order terms are reevaluated from the first global-pass 
solutions. These terms are held constant again during the solution for 
the second approximation obtained at the end of the second global pass. 
This procedure is repeated until the flowfield results corresponding to 
successive global passes converge within a specified limit [23]. 

4.5 Results and Discussions 

Numerical solutions of the viscous shock-layer (VSL) equations for 
the low-density hypersonic flow over long slender bodies are obtained. 
The surface slip [31] and the recently corrected shock-slip boundary 
conditions [62] are implemented in the implicit finite-difference method 
used to solve the governing equations. Detailed comparisons with the 
experimental data are included for several conditions. Extensive 
results are provided for long slender bodies with temperature surface 
conditions ranging from adiabatic to highly cooled. Also included are 
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Fig. 4.1 Solution sequence with shock and body slip. 
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calculations from the steady-state Navier-Stokes equations. These 
results provide an indication of the range of applicability of the 
viscous shock-layer solutions. 

4.5.1 Comparison with Experimental Data 

There are only a few experimental data for the low-density, high- 
energy flows available in the literature. The data of Little [46] are 
still considered quite good for such flows. These data, however, are 
limited to the measurements of pressure, drag, and skin friction. These 
data were used extensively for comparison with the theoretical 
predictions by Davis [12]. The same data have been employed for 
comparison with the present results also. Figures 4.2 to 4.6 give 

comparisons between the viscous shock-layer predictions and the 

0 

experimental data [46] for a 10 hyperboloid. With shock and body 

slip, surface pressure predictions by the present method agree quite 

well with the experimental data as shown in Fig. 4.2. Comparisons of 

0 

drag coefficients on a 10 hyperboloid for a range of values of the 
rarefaction parameter e are shown in Figs, 4.3 to 4.6. Predictions of 
Davis [12] are also given in these figures. It is clear that the 
present predictions with the shock and body slips are in much better 
agreement with the experimental data than the predictions of Davis. 

Large differences in the present calculations and those of Davis are 
seen with increasing values of e. These differences may be due to the 
errors in the slip conditions as mentioned in Sec. 4.1. 

Comparison between the predicted Stanton number distribution and 

0 

experimental data [64] for a 10 sphere-cone is provided in Fig. 4.7. 

The present calculations with shock and body slip are in good agreement 
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Fig. 4.3 Drag coefficient distribution for a 10* hyperboloid 
with e - 0. 1192. 
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with the data except for the stagnation point. The experimental heat 
transfer rate at the stagnation point were determined to be biased 
upward due to particle impact caused by the arc heater [65]. 

4.5.2 Comparison with Predicted Results 

Results obtained for a 22.5 hyperboloid in the stagnation region 
by the present viscous shock-layer equations and the steady-state 
Navier-Stokes equations are compared with those obtained by Anderson and 
Moss [23] in Table 4.1. The results from these calculations compare 
fairly well, especially for C p . A maximum difference of less than ten 

percent occurs between the two results at Rc,^ of 90 in the heat 
transfer coefficient, C H . This may be due to the grid clustering 

employed near the shock and body in the present calculations. 

Figures 4.8. and 4.9 show comparisons for the Stanton number and 
skin friction coefficient, respectively as obtained by the present 
viscous shock-layer method and that obtained by Davis [12]. Also shown 
are the present results obtained from the steady-state Navier-Stokes 
equations. The calculations are carried out for the stagnation point 
only for different values of the Reynolds number parameter, e, which is 
a measure of the degree of rarefaction. Larger values of e imply 
increased rarefaction effects. The two viscous shock-layer predictions 
have similar trends. However, significant differences are noticed for 
large values of e. The discrepancies may be due to the errors contained 
in the slip equations used by Davis [12] and the current grid clustering 
near the shock and body. Figures 4.8 and 4.9 also show that the viscous 
shock-layer predictions deviate from the Navier-Stokes results for large 
values of e. For e - 1, the present viscous shock layer predictions 
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give a Stanton number of 0.68, whereas the Navier-Stokes results yield a 
value of about 0.9. It is clear that the Navier-Stokes predictions have 
the right trend in approaching a value of 1 for the free-molecule flow 
limit. These results suggest that the viscous shock-layer approximation 

may not be valid for large values of e. 

Figures 4.10 and 4.11 show comparisons of predicted skin-friction 
coefficient and Stanton number distributions, respectively, as obtained 
by the present viscous shock-layer method and that by Gordon [16]. The 
comparison between these two results is quite good when the coarse grid 
structure of Gordon [ 1 6] is used. The method of Gordon is fully coupled 
and requires solving a 5 x 5 matrix at every point in the flowfield for 
a perfect gas. The complexity and stability problems in a fully coupled 
solution will be increased in analyzing a multi-species high-temperature 
air flow. Also, the computational times will be considerably large for 
long slender bodies by this method. The present approach, with coupling 
between the normal momentum and continuity equations only, may be more 
appealing for such flow conditions. Figures 4.10 and 4.11 also give 
results with and without slip for a variable grid near the shock and 
body surface. It is clearly seen that the computational grid-size as 
well as the slip effects are important in this case. 

11 , 5,3 Calculations for Different Altitudes and Surface Temperatures and 
Range of Validity of Viscous Shock-Layer Results 
An extensive test for the present computational method and the 
surface and shock-slip boundary conditions is provided through the 
results given in Figs. 4.12 and 4.13. The flow analyzed in these 
figures is a high Mach number (M^ - 20) flow over highly cooled (T w - 
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Flo. 4.9 Skin-friction coefficient at the stagnation point of 
22.5° hyperboloid as a function of epsilon. 
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Fifl. 4.10 Skin-friction coefficient distribution for a 22. 
hyperboloid with C - 0.223. 




Stanton number distribution for a 22.5* hyperboloid 
with e - 0.223. 
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300 K) , 2.54 cm nose radius, five and ten degree sphere cones. The 
free-stream conditions for very high-to-low altitude flight cases are 

given in Table 4.2. Figure 4.12 gives the Stanton number distribution 
0 

for a 5 sphere-cone obtained by using the viscous shock-layer equations 

for different body locations. The slip effects become insignificant at 

body locations greater than about ninety nose-radii or at altitudes less 

than about 60 km for a wall temperature of 300 K. The Stanton number 

0 

values with and without slip for a 10 sphere-cone as shown in Fig. 4.13 

0 

are higher than those for a 5 sphere-cone at the corresponding body 
locations for a given altitude except for the stagnation point (s - 0). 
At this location, the Stanton number values are almost the same for the 


O 0 

5 and 10 sphere-cones. Figures 4.12 and 4.13 also indicate that for a 
given altitude and body location, the slip effects are higher on the 


0 0 

conical flank portion for a 5 sphere-cone than for a 10 sphere-cone. 
Further, the effect of slip increases with increasing altitude for a 
given cone angle and body location. 

Figures 4.14 to 4.20 show the effect of surface temperature on 
stagnation-point pressure and heat-transfer coefficients. Both the 
viscous shock-layer and Navier-Stokes solutions are provided in these 
figures. The results presented in Figs. 4.14 to 4.17 show that the 
viscous shock -layer values of with no slip gradually increase from a 


value of 1.84 at about 30.5 km altitude to a value of 1.88 at 100 km 

altitude. The value of C stays constant at 1.84 for an adiabatic wall. 

P 


The viscous shock layer predictions for C 


P 


with slip continuously 


decrease with increasing altitude for a cooled surface. This trend is 
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300 K, V - 20.96 g/«ole, H - 20, S - 16.73 



Flo* 4.12 Stanton number distribution for a 5° sphere-cone as 
function of epsilon. 



Flo. 4.13 Stanton number distribution for a 10° sphere-cone as 
function of epsilon. 
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similar to the one given by the results of Davis [12]. This trend, 
however, indicates that viscous shock-layer results with slip do not 
approach the free molecule flow value at higher altitudes. The Navier- 
Stokes results with the body and shock slip do provide the right 
behavior of approaching the free molecule flow value at higher 
altitudes. These results, however, first show a decrease in the value 
of 0^ and then an increase as the altitude increases further. 

This behavior, also obtained by Jain and Adimurthy [66], is 
consistent with the trend observed by Potter and Bailey [67]. Good 
agreement between the Navier-Stokes results and the data of Potter and 
Bailey [67] was reported by Jain and Adimurthy [66], As can be noticed 
from Figs. 4,14 to 11.17, the dip in the pressure coefficient curve is 
reduced by increasing the wall temperature. For an adiabatic surface 
(implying no temperature slip), there is no dip in the C p curve, and it 

increases monotonically towards the free-moleeule flow value with the 
increasing in the altitude. It may be mentioned here that the free- 
molecule flow value of C p as well as its asymptotic value at lower 

altitudes is also influenced by the wall temperature. The free-molecule 
flow value is obtained from the equations of Bird’s [46]. The predicted 
value of C p from Navier-Stokes and viscous shock-layer solutions with 

slip approach the asymptotic value of 1.84, which is predicted by the 
inviscid modified Newtonian formulation, at lower altitudes with the 
increase in surface temperature. Obviously, this asymptotic value is 
obtained for a very high Reynolds number flow in absence of any slip 
effects. Reducing the slip effects by increasing the wall temperature 
also gives this asymptotic value at moderately high altitudes. 




Free 

molecule 
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Flo. *.15 Stagnation-point pressure coefficient for a 10 
sphere-cone versus Knudsen number and altitude 
with T,j - 1000 *K. 
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Flo. 4.16 Stagnation-point pressure coefficient for a 10 
sphere-cone versus Knudsen number and altitude 
with T* ■ 5000 *K . 
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Flo. 4.18 Staonatlon-polnt heat-transfer coefficient for 
a 10* sphere-cone versus Knudsen number and 
altitude with T* - 300 °K. 



119 



Fig. 4.19 stagnation-point heat-transfer coefficient for 
a 10* sphere-cone versus Knudsen number and 
altitude with T* - 1000 *K. 
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Fig. 4.20 stagnation-point heat- transfer coefficient for 
a 10* sphere-cone versus Knudsen number and 
altitude with T* - 5000 *K. 



121 


The results for the stagnation-point heat transfer coefficient as 
a function of the freestream Knudsen number (or altitude) are presented 
in Figs. 4.18 to 4.20 for different surface temperatures. Results for a 
given surface temperature show that the viscous shock-layer formulation 
with or without slip does not give physically realistic results at very 
high altitudes. The Navier-Stokes results with slip do approach the 
free-molecule flow value of approximately unity at very high altitudes. 
With an increase in surface temperature, the surface heat-transfer rate 
is decreased as expected. The effect of slip is noticeable down to an 
altitude of about 60 km for the various surface temperature considered 
here. Discrepancies of less than ten percent are noticeable below 
approximately 75 km altitude. 

Results of Figs. 4.14 to 4.20 suggest that the viscous shock -layer 
calculations with slip begin to deviate from the Navier-Stokes results 

ft ft 

with slip for freestream Knudsen number / R^) greater than about 

ft ft 

0.06. For 0^, the deviation begins at \ m / R^ * 0.01. It may not be 

appropriate to use the viscous shock-layer model even with body and 
shock slip at higher altitudes. 

4.6 Conclusions 

Results have been obtained for the surface pressure, drag, heat 
transfer, and skin-friction coefficients for hyperboloids and sphere- 
cone-shaped slender bodies under varying degrees of low-density flow 
conditions. Recently obtained surface-slip and corrected shock-slip 
conditions are employed to account for the low-density effects. The 
method of solution used for the viscous shock-layer equations is a 
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partially coupled spatial-marching implicit finite-difference technique. 
The flow cases analyzed include highly cooled surfaces in very high Mach 
number flows. The viscous shock-layer predictions compare quite 
favorably with experimental data. Results are also obtained from the 
steady-state Navier-Stokes equations by successive approximations by 
using the viscous shock-layer results to evaluate higher order terms for 
the first approximation. Comparison between the Navier-Stokes and 
viscous shock-layer results indicates that viscous shock-layer equations 
even with body and shock slip do not give physically consistent results 
in the stagnation region above approximately 75 km altitude for the 
conditions considered in this study. 
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Chapter 5 

CHEMICAL NONEQUILIBRIUM FLOWS 
5.1 Introduction 

One of dominant aspect of hypersonic flow is high-temperature 
effects. Strong compression of the gas forward of the vehicle and heat 
generation due to viscous dissipation lead to increase gas temperature 
in the shock layer. At high temperature, the gas will become chemically 
reacting. The specific heat per unit mass is considerably increased, 
the specific heat ratio will no longer equal to 1.4 and will no longer 
be a constant. The assumption of a calorically perfect gas is not 
appropriate; the effects of chemical reactions must be taken into 

account . 

From the example of atmospheric entry of the Apollo command 
vehicle given by Anderson [68], the shock layer temperature predicted on 
the basis of an equilibrium chemically reacting gas is a factor of 5 
less than the temperature predicted on the basis of a calorically 
perfect, nonreacting gas which gives an unrealistically high value of 
temperature in a high Mach number flow. Two major physical 
characteristics which cause a high-temperature gas to deviate from 
calorically perfect gas behavior, as stated in [68], are vibrational 

excitation and chemical reactions. 

All chemical processes take place by molecular collisions. As the 
temperature of the gas is increased, and hence the molecular collisions 
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become more violent, it is probable that the diatomic molecules of 
oxygen, 0 2 , and nitrogen, N 2> will be dissociated and nitric oxide, NO, 

will be formed or dissociated by collisions with other particles. In 
turn, collisions take time to occur. Hence, the chemical changes in a 
gas require a finite time to occur. Equilibrium flows assume that the 
gas has had enough time for the necessary collisions to occur. However, 
there are some flight conditions frequently encountered in atmospheric 
entry where the gas is not given the necessary time to come to a state 
of equilibrium. Under these conditions, flows are characterized by a 
chemical nonequilibrium process in the shock layer. The experimental 
wall temperature measurements and resulting heat-transfer rates obtained 
during the first flights of the Space Shuttle have been lower than the 
predicted equilibrium values at least over the first 40 percent of the 
Shuttle length and for much of the altitude range of interest [69-72] . 
The flight data from the Catalytic Surface Experiment [71,73] have 
verified that the lower rates can be attributed primarily to the fairly 
noncatalytic nature of the Shuttle thermal protection system and not to 
the unknowns in the freestream or flowfield quantities, and that the 
nonequilibrium effects persist to altitudes as low as 50 km for the 
orbiter . 

Although measurement data can be obtained from the Space Shuttle 
flights, it is very expensive for each flight. Moreover, a small scale 
laboratory experiment cannot simulate the chemical nonequilibrium flow 
around a hypersonic vehicle. An adequate design capability for future 
transportation systems relies on numerical predictions. Among the 
numerical methods available for solving the nonequilibrium flow over a 
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hypersonic vehicle, the viscous shock-layer approach not only has the 
advantage of requiring much less computing cost as compared to the 
Navier-Stokes solutions, but it also provides accurate predictions. 

This method has been used widely as a tool for engineering calculations 
[24,60,61 ,72,74,75] . 

The viscous shock-layer equations for a perfect gas and for a 
chemically reacting binary mixture were developed by Davis [12,56]. 

Based on this analysis, Moss [24] developed a code using the viscous- 
shock-layer equations for a multicomponent gas mixture with chemical 
equilibrium or nonequilibrium. It is shown that accurate results can be 
obtained by this nonequilibrium code [24,72,74,76,77]. However, the 
difficulties encountered in the case of a perfect gas are also 
encountered in the nonequilibrium flow over a slender long body. 

Appropriate shock and wall boundary conditions must be prescribed 
for the viscous shock-layer equations for chemically reacting flows. In 
addition to surface temperature and velocity, wall species 
concentrations are needed. However, the surface heating rate in a 
hypersonic nonequilibrium flow environment is strongly affected by the 
surface catalytic activity (the recombination of the dissociated oxygen 
and nitrogen atoms). For a dissociated nonequilibrium flow over a 
finite-catalytic wall, the heating rate and the wall species 
concentration is a function of the surface reaction rate coefficient (or 
energy— transfer recombination coefficient) [78]. The temperature- 
dependent and constant values of the coefficients for oxygen and 
nitrogen surface recombination have been determined by Scott [79]. 

These values have been incorporated in the viscous-shock-layer code 
available in [72]; however, the resulting heating predictions are only 
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in fair agreement with the STS-2 experimental data. Based on STS-2 
data, a new oxygen reaction rate expression has been developed by Zoby 
et al. [77], and it has been shown that a better heating comparison with 

experimental data can be obtained [75 ,773 • 

Thermodynamic properties and transport properties are required for 
each species considered in a multicomponent gas mixtures. All these 
properties are obtained from the polynomial curve-fit formulas based on 
measured data. Measurements made on the Orbiter during reentry have 
provided an extensive and reliable data base to improve these relations. 

The primary objective of this study is to investigate the effects 
of chemical nonequilibrium conditions in hypersonic flows over long 
slender bodies. For this, modifications in the existing code by Moss 
[24] are needed. The modifications included are:(1) the two first 
order equations, continuity and normal momentum, are solved 
simultaneously as a coupled set, and (2) the thermodynamic and transport 
curve fit relations are modified. The effects of different wall 
recombination coefficients on the predicted heat transfer are 
investigated. Also, this chapter includes a parametric study on the 
effects of body-angle, nose radius, and Mach number. 

5.2 Analysis 

The conservation equations that describe a reacting multicomponent 
gas mixture can be found in the literature [27,28]. The viscous shock- 
layer equations for nonequilibrium multicomponent gas mixture developed 
by Moss [24] are obtained from the conservation equations employing the 
same procedure as for the perfect gas [12]. For a blunt axisymmetric 
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body at zero angle of attack, the global continuity, s— momentum, and n 
momentum equations in the orthogonal, body-oriented transformed 
coordinates and in nondimensional form are written in the same form as 
Eqs. (2.25) through (2.29), (2.35) and (2.36) except the eddy viscosity, 

e + , is set to zero. For a chemical nonequilibrium flow, the energy 
equation is formulated in terms of the temperature instead of total 
enthalpy. In addition to these conservation equations, the species 
continuity equation and equation of state are needed to complete the set 
of equations. The energy and species continuity equations in the 
computational plane are written in the form of Eq. (2.25). The W in 
this equation represents T for the energy equation and for the 

species continuity equation. The coefficients to , Eqs. (2. M2) to 


(2.M9), are given here again by 
Energy, W - T: 

, . . n . k Jn . cose n N s 
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Species continuity, W - C.: 


1 3P W d 8 , "sh^ 

r- -35 ♦ h 


n .cose 
+ _ 3h 


°1 “ PL, 3n dTT M + nh 3h < r+ n h sh co30 


sh 


n ah n rti T<pu . n sh pv 
e 2 PL i (1+nn sh ic) e 2 PL 1 


(5.5) 


“ 2 “ 0 


(5.6) 


J_ 3PM i dg + ^i f "sh** 

nr 'N — . DI ' 


n .cose w.n 


sh 


) 


i sh 


3 " PL t 3n dn PL 1 M + nh sh < r+nn 3h cose^ e 2 pi ^ 


(5.7) 


n 3h pU 


£ PL. (1+nn . k) 
i sh 


where 


(5.8) 


PL i “ PF Ab ii 


(5.9) 



129 


PM 


i Pr 


N 3 

3C 

dg 

I S » lk 

K 

— 

k-1 

3n 

dn 

krfi 



Le, j 




( 5 . 10 ) 


i«k 


Ab 


ik 


(5.11 ) 


* N 

M i M ' 3 

Le >i - 1 ? u -ik * (l 
M 

jrfi 


M 3 

— ) I Le, C ] iltk 

K J- 1 J J 


In Eq. (5.11), Le * ik are the multicomponent Lewis numbers and M is 


molecular weight which is given by 
“* 1 
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The mass flux, , due to concentration gradients can be written as [29] 
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The equation of state is given by 

p - pT R / M* C p *. (5.14) 

The term w^ which appears in the energy and species continuity 
equations represents the rate of production of species i due to chemical 
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reactions. As discussed by Blottner [6,29] and Davis [56], the way of 
the production terms are written is very important in achieving 
convergence of the iteration procedure. Consequently, for the energy 
equation, the production terms are written as [56] 


( p ^k+1 


*. d w 

( "p + 3T ( ”p 5 k ( T k+1 ~ T k ) 


(5.15) 


where k denotes the iteration for which the solution is known and k+1 
the iteration for which a solution is required. It is found that if an 
expression of this type is not used which allows T to appear as an 
unknown in the energy equation, the method will not converge at low 
altitude conditions where the gas is approaching equilibrium conditions 
N 

[29]. The term £ s h.w, which appears in the energy equation, Eq. (5.3), 
i-1 1 1 

is written as follows: 

l* h i*i “ *1 + ™2 (5 ' 16) 


As for the species continuity equations, the production term is written 
as 


w 


-I - c i 4 i 


(5.17) 


Hence, Eqs. (5.15) and (5.17) express the production terms as a function 
of temperature for the energy equation and in terms of the species mass 
fraction for the species equations. 
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5.3 Boundary Conditions 

At the wall, the no-slip and no- temperature- jump boundary 
conditions are used in this study. No mass injection is considered, the 
normal component of velocity at the surface is taken to be zero. Also, 
the surface total enthalpy is given as 


N 

H - I 3 h C (5.18) 

X 4 ^ ^ 


For a nonequilibrium flow, the wall species concentration is 

* 

dictated by the catalytic recombination rate k y in the recombination 
equation which is given by 


aCj dg 
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(5.19) 


# £ 

where k , « k , / U . The catalytic recombination rate is determined 
w, i w, i ® 

from the catalytic recombination coefficient (or catalytic efficiency) 
T i by [78] 
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For a noncatalytic wall, the catalytic recombination rate of atoms is 
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equal to zero [313; hence, Eq. (5.19) becomes 


3C . 

‘ ST >w 


( 5 . 21 ) 


for all the species of a multicomponent mixture. For a fully catalytic 
wall, the catalytic recombination rate of each atom is equal to one 

[313. 

For a finite rate catalytic wall, the recombination coefficient 
for oxygen and nitrogen based on arcjet experimental heat-transfer data 

were determined by Scott [793 such as 
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(5.23) 


However, it has been noticed that the resulting heating predictions are 
only in fair agreement with the STS-2 experimental data by 
incorporating these coefficients [72,753. Zoby et al. [773 developed a 
recombination coefficient for oxygen based on experimental flight, STS- 

2, heat-transfer data. It is given by 

* 


, -658. 9/T 

Y Q - 0.00941 e 3 v 


(5.24) 


Boundary conditions immediately behind the shock are calculated by 
using the Rankine-Hugoniot relations. The nondimensional shock 
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relations are the same as Eqs. (3-21) to (3.23) and include 
Energy: 
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The chemistry across the shock wave is assumed to be frozen. 

5.4 Chemical Composition 

In this study, the chemical reaction is confined to a system of 
neutral air species (0, 0^, N, and NO). When chemical reactions 

proceed at a finite rate, the rate of production terms w^ are required. 

For a multicomponent gas with N g reacting chemical species and 

chemical reactions, the chemical equation describing the overall change 
from reactants to products may be written in the general form 
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where r - 1,2,..., N and N is equal to the sum of the reacting species 

r j 

(N ) plus the number of catalytic bodies. The quantities and B i>r 

3 

are the stoichiometric coefficients for reactants and products, 
respectively, whereas K* >p and K^ p are the forward and backward rate 

constants. The quantities X* denote the concentrations in moles per 

unit volume. The rate of change of any species as a result of a 
particular reaction is [80] 


* 

dX. 

<-i> 


- 7 > “ (B, r ~ <*,. r ) < K f r 11 X j 

dt r 1,r ’ ’ J-1 


J *°j.r 


N 


- K 


b ' P J-1 " J 


J *^1 r 

IT X . J ’ r ) 


(5.29) 


In order to find the net mass rate of production of the ith species per 
unit volume, Eq. (5.29) must be summed over all reactions r. Thus, the 

rate of production of chemical species, & lf can be expressed as 


* 


* 


i 


M, 


r 

I 

r«l 



dt r 


(5.30) 


The chemical reactions used in this study are as follows: 
r - 1 0 2 + M 1 - 20 + M 1 

N 2 + M 2 - 2N + M 2 


r - 2 
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r - 3 
r - 4 


N 2 + N - 2N + N 


NO + M 3 - N + 0 + M 3 


r - 5 
r - 6 


NO + 0 - 0 2 ♦ N 


N 2 + 0 - NO + N 


where M.^ , M 2 and M 3 are the catalytic third bodies [6]. The reaction 


constants for these equations are expressed in the modified Arrhenius 
form, where the forward rate is given as 


* # o 

« «C2 C x 10 J ' -a 

K f,r " T exp(iog e C0 r - JL— ) («i|) " 

T cm 


(4.31 ) 


and the backward rate is given as 

# * D2 r D1 x 10^ mole ~^r 

K b r ■ T “ ptl o* »° r - — ■ . ) .4- 

e T 3 cm 3 


(5.32) 


where 


j 

“r I “i r " 1 
i-1 ’ 


(5.33) 


and 


N. 

J 


l r - 1 

i-1 1,r 


(5.34) 


The values for the coefficients in Eqs. (5.31) and (5.32) are taken from 
the compilation of experimentally determined rate constants given by 
Blottner et al. [81]. For a specified temperature, density, and species 
composition, Eqs. (5.29) to (5.32) are used to determine the production 
rate of a multicomponent gas. 
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5.5 Thermodynamic and Transport Properties 
The thermodynamic properties C ^ and h^ and the transport 


properties and are required for each species considered. 


Since the multicomponent gas mixtures are considered to be mixtures of 
thermally perfect gases, the thermodynamic and transport properties for 
each species are calculated by using the local temperature and pressure. 
Then the mixture properties are determined in terms of the individual 
species properties. 

5 . 5.1 Thermodynamic Properties 

Values for the thermodynamic properties as a function of 
temperature are obtained by using polynomial curve fits for each 
chemical species. The following polynomial equations are used: 

Specific Heat 
# 



V 


a 3 T 


*2 


»o *4 

a,.T J ♦ a„T H 


(5.35) 


Enthalpy 


* * 
R T 


- a„ 


V 


a s T 


*2 


V 


V 


*4 


_6 

# 


(5.36) 


The development of these curve fits and a tabulation of the polynomial 
constant (a 1 to a^) are presented in [82]. These curve-fit formulas are 


tabulated up to 15,000 K. However, there are many flow conditions where 
the temperatures in the shock layer are much higher than 15,000 K. 

Hence, these curve fit formulas have been extended to a temperature 
range of up to 35,000 K by Shinn [83] based on the tabulated values 
given by Browne [84,85]* 
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5.5.2 Transport Properties 

Transport properties for viscosity and thermal conductivity are 
required for each species considered in the shock layer gas. These 
properties are obtained by using polynomial curve fits to the data of 
Esch et al. [86]. The mixture viscosity is obtained by using the 
semiempirical formula of Wilke [27]. The mixture thermal conductivity 
is obtained by a method analogous to that used for viscosity [24]. A 
binary diffusion model with Lewis number equal to 1.4 is used. 

In addition, the transport properties of the individual species 
are also obtained from the polynomial curve-fits in temperature to the 
values given by Yos [87]. These data are believed to be more accurate 
at the higher temperatures which are encountered in nonequilibrium 
calculations. With these individual species properties, transport 
properties for the gas mixture are obtained by using the methods of 
Armaly and Sutton [88] for viscosity and Mason and Saxena [89] for 
thermal conductivity. 


5.6 Method of Solution 

The method of solution is essentially the same as that used for 
solving the viscous shock-layer equations for one component perfect gas. 
The solution for the multicomponent gas mixture proceeds in exactly the 
same way as given in Chap. 3 for the one component gas except a species 
equation is included now. 

Three second-order equations, species continuity, s-momentum and 
energy, are replaced with central differences in the rrdirection and 
two-point backward differences in the ^-direction. The two first-order 
equations, continuity and n-momentum, are solved simultaneously. The 
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density in these two equations are replaced by Eq. (5.1 M). Then these 
two equations can be written in the same forms as Eqs. (3-60) to (3-63). 
The coefficients for a multicomponent gas mixture in these equations are 
given in Appendix A. All these equations along with the boundary 
conditions constitute a system of the tridiagonal form such as Eqs. 
(3.53), (3 . 64 ) and (3.65). The shock stand-off distance is evaluated by 
integrating the continuity equation and density is obtained from the 
equation of state. 

For specified free-stream conditions and body geometry, a 
stagnation streamline solution is obtained. With the stagnation 
streamline solution providing the initial conditions, the conditions at 
the shock providing the outer boundary conditions, and the conditions at 
the wall taken as the inner boundary conditions, the numerical solution 
is marched downstream to the desired body location At any body 
station m, the converged profiles at station m-1 are used as the initial 
guess for the profiles at station m. The solution is then iterated 
locally until convergence is achieved. The solution is advanced 
subsequently to the m+1 station. Figure 5.1 presents the procedure for 
solving the governing equations for any location m. 

5.7 Results and Discussions 

Numerical solution to the previously discussed viscous shock-layer 
equations with chemical nonequilibrium are presented and discussed. 
First, comparisons of present results are made with data of STS-2 [90] 
to investigate the effects of modifications of the chemical 
nonequilibrium code and the heating effects of different finite-rate 
oxygen surface recombination expressions. Second, results of three 
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different *U body-angle apnere cones are presented -bleb demonstrate 
the effects of surface catalysis and body angle. Finally, the effects 
of nose bluntness, Mach number and the thermodynamic and transport curve 
fit relations on surface heat transfer rate are Investigated. 

in the original code of Moss [24], the convergence criteria for 
each body station mas that the relative difference be less than 0.00, 
for both the temperature and tangential velocity derivatives at the 
wall. In this study, the temperature, tangential velocity, pressure and 
species concentration profiles are added into the convergence criteria 
where the relative difference is less than 0.00, for all these profiles. 
5.7.1 Comparison of the Present Method with Cascading Method 

Figures 5.2 and 5.3 show the results of shock standoff distance 

and surface heat transfer distribution over a 3b" sphere cone with 
finite rate chemistry. Without coupling the two first-order equations, 
continuity and normal-momentum, oscillation exists in the vicinity of 
the tangency point, this is the same as with the perfect gas model 
(Figs. 3.5 and 3.6). This oscillation can be removed by solving the two 

first-order equations simultaneously in a coupled way. 

5.7.2 comparison of the Present Method w ith Heasured Data 

Present predicted heating rates are compared with the STS-2 
laminar heating data [902 for three different altitudes in Figs. 5.4 to 
5.6. The viscous shock-layer equations are applied to the windward 
symmetry plane by using the concept of an equivalent axisymmetric bod, 
at zero degree angle of attack [72,77.91,922. Results for an altitude 
of 71.29 km are given in Fig. 5.4. With coupling the global continuity 
and normal momentum equations, the predicted heating rates are lower 
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Fig. 5.2 Shock stand-off distance for a 35* sphere-cone 
with nonequilibrium chemistry. 
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Fig. 5.3 Surface heat transfer rate distribution for a 
*35* sphere-cone with nonequilibrium chemistry. 
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than that without the coupling. Since the body angle is not small, the 
coupling effect is not significant. It is also shown that the predicted 
heating rates using the oxygen surface reaction rate expression of Zoby 
et al. [77] yield better comparison with the experimental data than that 
using the Scott’s relation [79]. The prediction using Scott's 
expression are 30 to MO percent lower than the experimental data. The 
present "equivalent” axisymmetric body results agree quite well with the 
three-dimensional results which are obtained by Thompson [75]. Figures 
5.5 and 5.6 give the comparisons of surface heating rate between 
predictions and STS-2 data at altitudes 60.56 km and 52.97 km, 
respectively. With finite catalytic wall conditions, the present 
predictions are in good agreement with the data over the length of the 
vehicle. 

The results presented in Figs. 5.M to 5.6 show that the flowfield 
in the shock layer at high altitude is quite far from that predicted by 
assuming the condition of chemical equilibrium. With decreasing 
altitudes, the data and the finite rate chemistry predict the results 
that approach the equilibrium value. However, near the stagnation 
region some degree of nonequilibrium flow persists to altitudes as low 
as 50 km, as shown in Fig. 5.6. 

5.7.3 Effects of Surface Catalysis and Body Angle 

O 0 

Results for three sphere cones with body half-angles of 20 , 10 , 

0 

and 6 are presented to illustrate the effects of body angle and surface 
catalysis. Freestream conditions are those for 53- 3^ km altitude and a 
Mach number of 25. The bodies have the same nose radius which equals 
0.0381 m. Both noncatalytic and fully catalytic surfaces are examined 
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Fig. 5.5 Comparison of predicted and experimental heating 
rate distribution at an altitude of 60.56 Km. 
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Comparison of predicted and experimental heating 
rate distribution at an altitude of 52.97 Km. 
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to show the limiting effects of wall catalycity on heating. Figure 5.7 

0 

gives the convergence history of the shock standoff distance for a 1 0 
sphere-cone using the present method. It does not matter whether the 
initial shock shape is created from the perfect gas solution for the 
same body-angle (the solid line and the dash line) or from the chemical 
nonequilibrium solution for larger body-angle (the chain line), the 

shock standoff distance will converge to the same shock shape. 

0 0 0 

The predicted heat-transfer distributions for 20 , 10 , and 6 
sphere cones are given in Figs. 5.8 through 5.13. In order to present 
these results clearly, it is necessary to show them in two figures for 
each body-angle. The first of these two figures is an enlarged view of 
the nose region, while the second extends up to 400 nose radii. From 
these figures, it is seen that the heating rates with the noncatalytic 
surface may decrease more than 50 percent in comparison to that with the 
fully catalytic surface on the spherical region. The differences in the 
heating rates decrease in the downstream regions. With a noncatalytic 
surface an appreciable amount of dissociation is present at the wall, 
then a diffusion-inhibiting blanket of unrecombined atoms can pile up 
near the surface and thereby reduce the heat transfer to the wall by 
diffusion; this could result in a reduction of heat transfer to the 
surface. On the other hand, with a fully catalytic surface the atomic 
species which result from the dissociation in the high temperature air 
will recombine with oxygen and nitrogen molecules at the surface. Since 
the diffusion flux of atoms toward the wall must be equal to the rate of 
disappearance of atoms at the wall due to recombination there, the heat 
transfer to wall by diffusion is maximum with a fully catalytic wall. 
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Fig. 5.7 Shock stand-off distance for a 10* sphere-cone 
with nonequi 1 ibrium chemistry. 
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Flo. 5.8 Nonequilibrium surface heat transfer rate 

distribution for a 20“ sphere-cone (fore-cone). 
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Fig. 5.9 Nonequilibrium surface heat transfer rate 

distribution for a 20 * sphere-cone (aft-cone). 



Alt=53.34 km 
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Flo. 5.10 Nonequl 1 ibrium surface heat transfer rate 

distribution for a 10* sphere-cone (fore-cone). 
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Fig. 5.11 Nonequi librium surface heat transfer rate 

distribution for a 10° sphere-cone (aft-cone). 



Alt=53.34 km 
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Fi g . 5.12 Nonequi 1 ibrium surface heat transfer rate 

distribution for a 6 * sphere-cone (fore-cone). 
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For a real surface, some degree of catalysis is usually present, hence, 
the heating rate for a finite catalytic surface is between these two 
limits. 

The results for the ratio of surface heating rate with 
noncatalytic wall to that with fully catalytic wall for the three bodies 
are shown in Figs. 5.14 and 5.15. The ratio demonstrates the maximum 
potential for a surface heating-rate reduction in the presence of 
dissociated nonequilibrium flow over a finite-catalytic surface. All 
the three curves have the same trend. The ratios keep decreasing up to 
the tangency point, then increase up to a maximum value in the 
recompression region, and finally, decrease to a constant value on the 
far downstream region. It is noticed that the location of the maximum 
point moves downstream as the body-angle decreases; i.e. , the 
recompression region moves downstream as the body-angle decreases. 

For stations beyond 100 nose radii (Fig. 5.15), the results 

0 

indicate less nonequilibrium effects for the slender 6 cone than the 
0 0 

10 and 20 cones. As shown in Figs. 5.16 and 5.1.7, these can be 
attributed to more dissociated species present throughout the flowfield 
for energy transport by diffusion to the surface for the wider angle 
cone. In this body region for the lower cone angles, conditions 
sufficient to produce dissociated species exist only in a small region 
of the boundary layer. In the fore-cone region (Fig. 5.14), the largest 
cone angle produces the smallest nonequilibrium effects as indicated by 
higher values of the ratio. (This result does not imply the local or 
total heating rate to the larger cone angle is less.) As shown in Figs. 
5.18 and 5.19, this trend may be explained for a given nose radius, the 
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flow over the smaller cone angle compared to the larger cone angle 
expands rapidly and results in freezing of the flow chemistry and larger 

percentages of dissociated species. 

Figures 5.20 to 5.22 show the effects of body-angle on the surface 
heating rates and wall pressure. Figures 5.20 and 5.21 give the heating 
rates for a noncatalytic wall and for a fully catalytic wall, 
respectively. The wall pressure distributions are illustrated in Fig. 
5.22. Since the wall catalyticity has negligible effect on the wall 
pressure, the results are presented only for the noncatalytic wall. 
Decreasing the body angle can reduce the wall pressure and heat transfer 
to the surface. Hence, the desirable geometries for hypersonic vehicles 
are slender long bodies in order to reduce the heat transfer rate and 

drag force on the bodies. It is noticed in Fig. 5.22 that the 

* * 0 
recompression regions start at s /I^ - 3.0, 10.0, and 28.0 for the 20 , 

10° and 6° sphere cone, respectively, under the condition investigated. 
5 . 7 .H Effects of Nose Bluntness and Mach Number 

A study of the effects of nose blunt ness and Mach number on the 
shock standoff distance, surface heat transfer and mass concentration of 

0^ at the stagnation point of a 45 sphere cone is conducted and results 

are discussed here. Freestream conditions are selected at 90 km, nose 
radii are 0.305 m to 2.286 m, and Mach numbers are 30 to 36. Only 
noncatalytic surface boundary condition is considered. 

Figures 5.23 to 5.24 show the effects of nose blunt ness on the 
shock standoff distance and surface heating rate for a Mach number of 
36 . It is seen that the shock standoff distance increases and surface 
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F1 °* 5.15 Effect of body angle on noncatalytlc to fully 

catalytic heating ratio distribution (aft-cone). 
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Fig. 5.16 Species concentration profiles at s - 290.0 for a 
noncatalytic wall. 



Fig. 5.17 Species concentration profiles at s • 290.0 for a 
fully catalytic wall. 
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heating rate decreases as the nose radius increases. Figure 5.25 shows 
the concentration profile of 0 2 for three different nose radii. Since 

the shock standoff distance increases as the nose radius increases, both 
the oxygen and nitrogen molecules have more time to dissociate. 
Consequently, the atomic oxygen and nitrogen concentrations are higher 
in the shock layer. Moreover, the amount of energy absorbed by the 
dissociation phenomena has reduced the temperature in the shock layer 
such that the heat transfer to the wall is decreased. 

Figures 5.26 and 5.27 show the effects of Mach number on the shock 
standoff distance and surface heating rate for a nose radius equal to 
0-91^ m. The shock standoff distance decreases and the surface heating 
rate increases as Mach number increases. At higher Mach numbers, the 
density and pressure increase across the shock wave are larger and hence 
the mass flow behind the shock wave can readily squeeze through a 
smaller area. Moreover, the temperature increase across the shock wave 
is larger as Mach number increases. As a result, a larger amount of 
energy transfers in the shock layer; i.e., the amount of heat transfer 
to the surface is larger. Figure 5.28 shows the concentration profile 
of oxygen molecules for three different Mach numbers. At the higher 
Mach numbers, oxygen dissociation increases due to the larger 3hock 
layer temperature. 

It is noted from Fig. 5.26 that the shock standoff distance 
changes less and less as Mach number increases; i.e., the shock standoff 
distance becomes relatively insensitive to changes in free stream Mach 
number at high Mach numbers. This is an example of the Mach number 
independence principle. 
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Effect of body angle on nonequl librium heating 
rate with noncatalytic wall. 


Fig. 5.20 
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Fig. 5.21 Effect of body angle on nonequilibrium heating 
rate with fully catalytic wall. 
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Flo. 5.22 Effect of body 3nsl6 on surfsce pressure 
distribution. 
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Fig. 5.23 Effect of nose radius on shock stand-off distance 
at the stagnation point. 
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Fig. 5.25 Effect of nose radius on concentration 

distribution of O at the stagnation point. 
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Fig. 5.26 Effect of Mach number on shock stand-off distance 
at the stagnation point. 
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5.7.5 Effects of the Thermodynamic and Transport Curve Fit Relations 

The polynomial curve-fit formulas for thermodynamic and transport 
properties in [72] are based on Esch’s data [82,86]. The mixture 
viscosity and thermal conductivity is obtained by Wilke’s formula [27]. 
Values for the coefficients in these formulas are tabulated only up to 
15,000 K. However, under certain conditions the temperatures in the 
shock layer are considerable higher than 15,000 K. Therefore, new 
polynomial curve-fit formulas for species transport properties have been 
developed from the data of Yos [87,93]. The gas mixture transport 
properties are obtained using the method of Armaly and Sutton [88] for 
the viscosity and Mason and Saxena [89] for the thermal conductivity. 

The polynomial curve-fit formulas for thermodynamic properties have been 
extended to an upper temperature range up to 35,000 K by Shinn [83]. 
These provide accurate predictions at higher temperatures encountered in 
nonequilibrium calculations [26, 74,75]. The Prandtl number in the shock 
layer is set equal to a constant [24,26] or is calculated from the local 
thermodynamic and transport properties [75,76]. A study of the effects 

of the thermodynamic and transport properties on the surface heating 
0 

rates of a 6 sphere cone is presented. Freestream conditions are those 
for an altitude of 53-34 km, a nose radius is 0.0381 m and a Mach number 
is 25. A noncatalytic surface boundary condition is considered. 

Surface heating rate results with Yos' as well as with Esch's 
transport property data are given in Fig. 5.29. The results obtained by 
Thompson [34] are also shown in the figure for comparison. The Prandtl 
number is computed using local thermodynamic and transport properties. 




Flo. 5.29 Effect of the transport curve fit relations on surface 
heat transfer rate for a 6" sphere-cone. 



Alt - 53.34 Km 



Flo. 5.30 Effect of the transport curve fit relations on surface 
heat transfer rate for a 6* sphere-cone (fore-cone). 
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Effect of the transport curve fit relations on surface 
heat transfer rate for a 6* sphere-cone (aft-cone). 
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The maximum differences between the two sets of results are about five 
percent. The differences between Thompson's results and the present 
results, which are obtained using Yos* data, are about nine percent. 

The reason for the discrepancy may be due to using different shock shape 
and normal-direction grid size. The outer boundary conditions depend on 
the shock shape, and the gradients of the flowfield quantities depend on 
the grid size. 

Figures 5.30 and 5.31 show the surface heat transfer rates for 
variable as well as constant Prandtl number. Both Yos' and Esch's data 
have been used to compute the individual transport properties. The 
constant Prandtl number is set equal to 0.72. The predictions for the 
constant Prandtl number with Esch's data are always lower than the other 
predictions. On the other hand, the predictions for the variable 
Prandtl number with Yos' data are the highest values. Using Yos' data, 
higher surface heat transfer predictions are obtained for constant or 
variable Prandtl number than that using Esch's data. The differences 
between that using the variable Prandtl number with Yos' data and that 
using a constant Prandtl number with the Esch's data are about ten 
percent in the recompression region and about five percent on the other 
locations . 


5.8 Conclusions 

Numerical solutions of the viscous shock-layer equations under 
nonequilibrium chemistry conditions are presented for hypersonic flow 
over long slender bodies. The method of solution used for the viscous 
shock-layer equations is a partially coupled spatial-marching implicit 
finite-difference technique. The flow cases analyzed include 
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noncatalytic, finite catalytic, and fully catalytic surfaces. Results 
from the present method show that the coupling between the global 
continuity and normal momentum equations is essential and adequate to 
obtain stable solutions past long slender bodies. The comparisons of 
the present predictions with the STS-2 laminar heating data indicate 
that the oxygen surface reaction rate expression of Zoby can give better 
agreement with the flight data than the extrapolation of ground based 
experimental recombination data. It is shown that near the stagnation 
region of the vehicle, some degree of nonequilibrium flow persists to 
altitudes as low as 50 km. 

It is shown that surface catalytic effects as well as body angles 
Influence significantly the surface heat transfer rates. For a 
noncatalytic surface, the heating rates can decrease more than 50 
percent in comparison to that for a fully catalytic surface near the 
stagnation region. These effects become less significant in the 
downstream region. It is also shown that the heating rate due to 
diffusion for a smaller body-angle sphere cone is not as important as 
for a larger body-angle in the downstream region. 

In order to reduce surface heat transfer rate and drag force on 
the bodies, the desirable geometries of hypersonic vehicles are slender 
bodies with blunt noses. Although increasing the nose bluntness can 
decrease the surface heating rate, it will increase the pressure drag 
coefficients. Optimization should be made between the drag and heat 
transfer rate. 

With thermodynamic and transport properties from Esch's data and 
for a constant Prandtl number, the present method always predicts lower 
heating rate than that from Yos' data. For higher temperature 
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conditions, the polynomial curve-fit formulas for species transport and 
thermodynamic properties based on Yos* and Browne's data can give better 
predictions of surface heat transfer. 
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Chapter 6 

CONCLUDING REMARKS 

A method for solving the viscous shock-layer equations for 
hypersonic flows over long slender bodies is presented. These equations 
are solved by employing a spatial-marching, implicit finite-difference 
technique. The two first-order equations, continuity and normal 
momentum, are solved simultaneously as a coupled set. This method 
yields a simple and computationally efficient technique. 

A wide range of flow conditions has been considered in this study. 
This includes conditions from high Reynolds number at low altitudes to 
low Reynolds number at high altitudes. At low altitudes, the hypersonic 
flow over a slender body usually becomes turbulence. Two algebraic 
turbulence models, Cebeci-Smith and Baldwin-Lomax, have been used with 
the present numerical technique for application to long slender bodies. 
At high altitudes, the low density effects become important. Recently 
obtained surface-slip and corrected shock-slip conditions are employed 
to account for these effects. At higher altitudes, the gas becomes 
chemically reacting. Under certain conditions, the flows are 
characterized by chemical nonequilibrium conditions in the shock layer. 
Numerical solutions under these conditions are also obtained for long 
slender bodies. 

Results for different conditions are obtained for axisymmetric 
bodies at an angle of attack of zero degree. Detailed comparisons have 
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been made with other predictions and experimental data for slender body 
flows to assess the accuracy of the present numerical technique. The 
results show that the coupling between the continuity and normal 
momentum equations is essential and adequate to obtain stable and 
accurate solutions past long slender bodies. This is true for both the 
chemically nonreacting and reacting flows. 

It is shown that both the Cebeci-Smith and Baldwin-Lomax models 

are adequate for application to long slender bodies. Due to the 

sensitivity of the Cebeci-Smith turbulence model to the boundary layer 

edge location, it is imperative that the numerical method provide good 

resolution and accurate solutions near the boundary layer edge. The 

Baldwin-Lomax turbulence model, which avoids the use of conventional 

boundary layer thickness in its formulation, appears more convenient to 

implement. However, it has been shown that the constant, C , depends 

cp 


on the flow Mach number. Based upon this study and other 
investigations, a linear dependence of this constant with Mach number is 
suggested. Further comparisons with experimental data are needed to 
verify this dependence. 

Using the corrected slip models, the viscous shock-layer 
predictions compare quite favorably with experimental data. The slip 
effects become insignificant in the downstream region or at altitudes 
less than about 60 km for geometry and conditions considered in this 
study. Significant slip effects are observed primarily in the 
stagnation region. Comparison between Navier-Stokes and viscous shock- 
layer results indicates that viscous shock-layer equations, even with 
body and shock slip, do not give physically consistent results in the 
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stagnation region above approximately 75 km altitude for the conditions 
considered here. 

The present prediction with finite-rate chemistry yields good 
comparison with the STS-2 laminar heating data. Near the stagnation 
region of the vehicle, some degree of nonequilibrium flow can persist to 
altitudes as low as 50 km. Under chemical nonequilibrium conditions, 
the surface catalytic effects can influence significantly the surface 
heat transfer. For a noncatalytic surface, the heating rates can 
decrease more than 50 percent in comparison to that for a fully 
catalytic surface near the stagnation region. These effects become less 
significant in the downstream region. The heating rate due to diffusion 
for a smaller body-angle sphere cone is not as important as for a larger 
body— angle sphere cone in the downstream region. 

For further study, it is recommended that the present method be 

used to study the following physical problems 

1 . incorporate surface and shock slip conditions in the finite 
rate chemistry code to investigate the low density effects. 

2. increase the number of species, including ionized species. 

3. incorporate a convenient radiative transfer model to 
investigate the effects of radiative heat transfer in the 
shock layer. 

4. analyze the effects of angle of attack by developing three- 
dimensional viscous shock-layer equations. 

5. modify the viscous shock-layer equations with equilibrium 
chemistry for long slender bodies. 
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APPENDIX A 


COEFFICIENTS OF CONTINUITY AND NORMAL MOMENTUM 
EQUATIONS AS A COUPLED SET 


The coefficients of Eqs. (3.60)-(3.65) are given in this appendix. 
For a perfect gas, the coefficients in Eqs. (3.60)-(3.63) are given as 
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The coefficients K 1 to K g in Eqs. (3-64) and (3-65) are given by 


K 1 “ (A c,n+1/2 C NM, n+1/2 * A NM,n+1/2°c t n+1/2 ) (B c,n-1/2 A NM,n-1/2 


B NH,n-1/2 A c,n+1/2 ) 


(A. 14) 


K 2 “ (A c,n+1/2°NM, n+1/2 " A NM,n+1/2 D c,n+1/2 ) (B c,n-1/2 A NM, n-1/2 


B NM,n-1/2 A c,n-1/2* “ ^ A c, n+l/B^M, n+1/2 ” A NM,n+1/2 B c,n+1/2^ 



196 


x (B c,n-1/2 C NM, n-1/2 ” B NM,n-1 /2 C c,n-1 /2^ 

K 3 “ " (A c,n+1/2 B NM,n+1/2 " A NM,n+1/2 B c,n+1/2 )(B c,n-1/2 D NM 
B NM,n-1/2 D c,n-1/2 ) 

K 4 (A c, n+ 1 / 2 E NM , n+ 1 / 2 ~ A NM,n+1/2 E c,n+1/2 )(B c,n-1/2 A NM 

B NM,n-1/2 A o, n-1/2* " (A c,n+1/2 B NM,n+1/2 “ A NM,n+1/2 B c 
X (B c, n-1/2 E NM,n-1/ 2 ~ B NM,n-1/2 E c, n-1/2* 

K 5 " (C c,n+1/2 A NM,n+1/2 ~ Sw,n+1/2 A c,n+1/2* (D c,n-1/2 C NM,n- 
“ °NM,n-1/2 C c, n-1/2* 

K 6 " (C c,n+1/2 B NM,n+1/2 " Sw,n+1/2 B c,n+1/2* (D c,n-1/2 C NM,n 
D NM,n-1/2 C c, n-1/2* ~ (C c,n+1/2 D NM,n+1/2 " °NM,n+1/2 D c 
X (D c,n-1/2 A NM,n-1/2 ” D NM,n-1/2 A c, n-1/2* 

K 7 " " (C c,n+1/2 D NM,n+1/2 “ C NM,n+1/2 D c,n+1/2 )(D c,n-1/2 B NM 


(A. 15) 

n-1/2 
(A. 16) 

n-1/2 " 

n+1/2* 
(A. 17) 

1/2 

(A. 18) 

1/2 “ 

n+1/2) 
(A. 19) 

n-1/2 



197 


D NM,n-1/2 B c,n-1/2 ) 


(A. 20) 
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APPENDIX B 


NAVIER-STOKES EQUATIONS IN THE BODY ORIENTED 
COORDINATE SYSTEM 


The steady-state form of the Navier-Stokes equations is taken from 
Anderson and Moss [233. For an axisymmetrlc or two-dimensional body at 
zero angle of attack, these equations In the body-oriented coordinate 
can be written in the form of Eq. (2.25). All the coefficients in Eqs. 
(2.26) to (2.34) are the same except Eqs. (2.28), (2.32) and (2.36) are 
need to modify. This modification is given in this appendix. 

For the steady-state Navier-Stokes equations, coefficient 

appearing in Eq. (2.25) should be replaced by a' 3 , which is defined as: 
s-momentumu 
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where a 3 in Eqs. (B.1) and (B.2) is given by Eqs. (2.28) and (2.32), 

respectively. Abbreviation HOT in Eqs. (B.1) and (B.2) represents 
higher-order terms, and subscripts 's' and 'e' imply terms in the s- 
momentum and energy equations, respectively. 

The normal momentum equation, Eq. (2.36), for the Navier-Stokes 
model is modified to 

<LHS V (2.12) - (H0T) n- 0 (B - 3> 

where the first term in Eq. (B . 3) implies the entire left-hand side of 
Eq. (2.36), and the second term represents the higher-order terms in the 
n-momentum equation. 

The higher-order terms appearing in Eqs. (B.1) through (B.3) are 
defined as: 
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Variables through C appearing in Eqs. (B.H) through (B.6) are given 


by : 
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